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NAVIER-STOKES ANALYSIS OF AIRFOILS WITH LEADING EDGE ICE ACCRETIONS 

Mark G. Potapczuk 
University of Akron 
Akron, Ohio 44325 

ABSTRACT 

A numerical analysis of the flowfield characteristics and the performance dograda 
of an airfoil with leading edge ice accretions was performed. The important fluid dyn; 
processes were identified and calculated. Among these were the leading edge separation 1 > u 
at low angles of attack, complete separation on the low pressure surface resulting in prema 
stall, drag rise due to the ice shape, and the effects of angle of attack on (lie separated 
field. Comparisons to experimental results were conducted to confirm these calculations. 

A computer code which solves the Navier-Stokcs equations in two dimensi 
ARC2D, was used to perform the calculations. A Modified Mixing Length turbulence 111 
was developed to improve capabilities in calculating the separated flow phenomena. A 
generation code, GRAPE, was used to produce grids for several ice shape and ai 
combinations. 

Results indicate that the ability to predict overall performance characteris'd irs. 
as lift and drag, at low angles of attack is excellent. Transition location is important 
accurately determining separation bubble shape. Details of the flowfield in and downstreai 
the separated regions requires some modifications. Calculations for the stalled airfoil indi 
periodic shedding of vorticity that was generated aft of the ice accretion. Time aver, 
pressure values produce results which compare favorably with experimental informal. ion 
turbulence model which accounts for the history effects in the flow may bo justified. 
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CHAPTER 1 


INTRODUCTION 

This work details the investigation of the aerodynamics of an airfoil with leading edge 
ice accretions. A computational approach is employed for the prediction of the fluid dynamics 
and performance characteristics of such an airfoil. The method presented is based on the 
numerical solution of the Navier-Stokes equations in a body-fitted coordinate system. A 
previously existing airfoil code, along with a new turbulence model developed for this study, 
was adapted to incorporate the physics of a leading edge separation bubble and to account for 
the completely separated flow at high angles of attack. This code was used to examine the 
structure of the separation bubble, the development of the turbulent boundary layer aft of the 
bubble, and examine the effects of angle of attack on the stalled airfoil flowfield. In addition, 
the code was used to determine the changes to maximum lift, angle of attack at stall, and drag 
as a result of the presence of ice on the leading edge. 

Comparisons to experimental results, for several airfoils with artificial leading edge 
ice shapes, are presented to verify the method and to illustrate its predictive capability. 
Capabilities for the quantitative analysis of airfoil performance degradation due to icing were 
previously limited to experimental correlations and potential flow analysis of relatively 
aerodynamic ice accretions. Evaluation of less benign ice shapes requires the ability to calculate 
separated flow regions and the resulting changes in lift, drag, and moment forces. This can be 
accomplished by the use of either an interactive boundary- layer approach or the solution of the 
Navier-Stokes equations. The former method is under evaluation by Cebeci [1]. The 
investigation of the use of a Navier-Stokes code thus seems justified in order to examine the 
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relative merits of each approach. Optimum design of ice protection systems generates the 
need for such computational capabilities. If the performance of an iced airfoil can be 
accurately predicted, then more meaningful design parameters can be considered in the 
development of ice protection systems. Further, if the aerodynamics of the iced airfoil are 
known then other analysis methods can be used to predict the ice buildup on the airfoil. This 
should serve to decrease the need for actual wind tunnel testing. 

1.1 Description of Problem 

Aircraft icing occurs when an aircraft passes through a cloud of supercooled water 
droplets and droplet impingement combined with heat transfer processes result in accretion of 
an ice mass which may cause significant performance loss to the aircraft. The severity of the 
icing encounter is dependant on meteorological conditions, flight conditions, aircraft geometry, 
pilot performance, and anti/de-icing measures employed. Prediction of aircraft capabilities 
under icing conditions is determined by the first three of these factors. 

Meteorological conditions determine the charcterisitics of the ice itself. Typically the 
types of ice identified for these purposes are rime, glaze, and mixed. Temperature is the major 
influence on determining the type of icing encountered. At temperatures just below freezing, a 
clear granular ice, termed glaze ice, tends to form which can produce very non- aerodynamic 
shapes on the airfoil leading edge. At very cold temperatures, an opaque relatively smooth ice, 
termed rime ice, forms which has a more aerodynamic shape than glaze. At some intermediate 
temperatures, a mixture of these conditions occurs usually with a glaze ice core surrounded by 
a shell of rime. These ice types combined with flight conditions, length of icing encounter, 
and aircraft geometry determine the final ice shape and resulting aerodynamic degradation 
characteristics. 
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The presence of ice on the leading edge surfaces of an aircraft can lead to severe 
degradation of the aircraft performance as a result of several influences. The weight of the ice 
itself can require additional lift to maintain the desired altitude. At the same time, the ability 
of the lifting surfaces to provide that lift is reduced. Additionally, the ice accretions increase 
the drag forces on the aircraft resulting in increased thrust requirements. Any one of these 
influences can be enough to prohibit operation of an aircraft in an icing environment. 
Evaluation of the performance degradation due to these influences is important for specification 
of appropriate ice protection measures. 

The weight of the ice can be accounted for as an additional payload. Thus, it should 
be within the scope of present design techniques to determine this influence. The increased 
drag and decreased lift due to the ice accretion are not so easily determined. Drag calculations 
require an ability to evaluate both pressure drag and skin friction. This means being able to 
determine the appropriate pressure distribution on a highly irregular geometry and to 
accurately determine the influence of the turbulent boundary layer on that surface. Lift 
calculations also require the ability to accurately determine pressure on the surface. These 
capabilities are available in existing codes only for attached flows on normal airfoil shapes. 

Airfoils with leading edge ice shapes pose greater difficulties for computation due to 
the flow characteristics associated with the non-aerodynamic geometry of the ice accretion. 
The ice shape produces a leading edge separation bubble which at low angles of attack 
reattaches to the surface. The size of this region and the flowfield inside it are difficult to 
model due to the limitations of the method of modeling turbulence previuosly employed in 
most computer codes. The location of the reattachment point is also affected by the 
turbulence model, which in turn can influence the development of the reattached boundary 
layer downstream of the bubble. At high angles of attack, this bubble detaches from the 
surface and results in an unsteady flowfield characterized by periodic vortex shedding. This 
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leading edge separation results in premature stall of the airfoil. Accurate modelling of these 
phenomena is essential for determination of the performance characteristics of the airfoil during 
an icing encounter. 

1.2 Icing Analysis Plan 

The development of an iced airfoil analysis code is part of a national icing analysis 
effort coordinated by NASA and the FA A. In this effort, the goal is to be able to calculate the 
growth of ice on an airfoil and determine the aerodynamic performance degradation due to 
that ice. Additionally, it is desired to be able to model various de-icing and anti-icing systems 
in order to limit expensive icing tunnel tests. This approach requires the development of 
several types of codes and ultimately the tying together of these codes in an overall analysis 
scheme. The codes used for icing analysis are employed in a computational loop which starts 
with evaluation of the flowfield for a clean airfoil, moves to a calculation of the water droplet 
trajectories, calculates the ice buildup from those particles which impinge on the surface, and 
evaluates the performance degradation due to this ice growth. The new shape is then used to 
re-evaluate the flowfield and go through the loop again. This process would be repeated until 
the ice buildup reached some critical level, most likely the stall point. The computational 
approach developed in this study would be used for calculation of the flowfield at the start of 
the loop and for performance evaluation at the end of the loop. 

For the overall icing analysis plan, it is essential, to be able to perform the 
calculations of the type reported in this work. The methods employed up to the present have 
either been unsuccessful or successful over a limited range of conditions. The use of a Navicr- 
Stokes code with an appropriate turbulence model holds the promise of achieving this critical 
capability for icing analysis research. 
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CHAPTER 2 


ICED AIRFOIL PERFORMANCE ANALYSIS : STATE OF THE ART 

The ability to predict the performance degradation of an airfoil due to ice deposition 
is not presently available to the aerospace community. Several empirical studies have been 
completed over the years, but no well-developed method has been produced to determine the 
loss in lift and increase in drag associated with this phenomena. The advent of high speed 
computers and development of efficient numerical schemes for solving partial differential 
equations has provided the opportunity for further development of prediction methods for iced 
airfoil aerodynamic analysis. It is the objective of this investigation to utilize these recent 
advances in computational ability and to examine their strengths and weaknesses in relation to 
the icing problem. This chapter presents a brief examination of the state of the art in icing 
analysis up to the present. The development of computational methods in the codes employed 
during this investigation is also discussed. 

2.1 Historical Background 

The current method of ice protection, hot air bleed, was first developed in the 1940’s 
as a result of an experimental program conducted by NACA [2]. This program consisted of 
numerous wind tunnel tests of iced airfoils, in an attempt to develop correlations for the 
associated drag rise. The correlations of Gray [3] were developed during this period and have 
served as the basis for predicting drag rise up until the present. Gray’s correlation is known to 
be a very approximate method and is considered inadequate for the design of ice protection of 
modern general aviation aircraft and rotorcraft. 
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During the past five to ten years, work in the area of correlation development has 
been re-examined in an attempt to improve upon Gray’s results. In 1982, Bragg [4] developed 
a method which accounted separately for the effects of changes in leading edge geometry and of 
surface roughness for rime ice profiles. He states that this method is limited to small ice 
accretions at low angles of attack. In 1983, another empirically based approach was attempted 
by Miller, Korkan, and Shaw [5]. Their conclusion was that further work was required for the 
development of a more general drag rise correlation. Additionally, these methods are used 
strictly for the prediction of global characteristics of the iced airfoil flow field. Detailed 
evaluation of the velocities, pressures, temperatures, or other important parameters, is not 
provided by the use of correlations. However, this type of information is necessary input to 
any method designed to predict ice accretion and growth, as envisioned for a comprehensive 
icing analysis program. 

Analytical methods for predicting drag rise due to icing were evaluated by Peterson 
and Dadone [6] in 1980. Their conclusion was that then current methods underpredicted drag 
rise. Later, in 1982, Bragg, Gregorek, and Shaw [7] used the ai rfoii code of Epplcr [8] to 
successfully predict the maximum lift coefficient mar ) f° r a rime ice profile. In 1984, 

Bragg [9] investigated the aerodynamic characteristics of airfoils with rime and glaze ice 
accretions. He compared several computer codes, eventually settling on the Dvorak CLMAX 
code [10], II is results for rime ice were encouraging, but several difficulties were found in 
calculation of drag and in prediction of leading edge separation on glaze ice profiles. Bragg 
also used the Bristow potential flow code [11] and measured separation bubble geometry to 
predict Cp values on a glaze ice shape. This approach does not, of course, allow evaluation of 
viscous effects. 

As mentioned previously, either an interactive boundary layer method or solution of 
the Navier-Stokes equations is necessary to produce the desired results. The interactive 
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boundary layer approach currently being used by Cebeci [l], shows promising preliminary 
results, especially for angles of attack below stall. However, some details of the flowfield, such 
as the reverse flow region of the separation bubble, are not predicted properly. Further, the 
massive separation which occurs at and above stall angles is not capable of being predicted by 
this method. Hence, it seems justified to employ an alternate approach, in an effort to 
complement the interactive boundary layer method, for situations which require the additional 
physical modeling provided by the Navier-Stokes equations. Thus, it may be envisioned that 
the less computationally expensive interactive boundary layer approach would be employed for 
small angles of attack and the Navier-Stokes solver would be used for the evaluation of stall. 
Additionally, the Navier-Stokes solver could be employed for evaluation of velocities within the 
separation bubble region. It was with these goals in mind that this work was undertaken. 

2.2 Background of Numerical Procedures 

The literature on grid generation, Navier-Stokes analysis, and turbulence modeling is 
quite extensive. Thus, any presentation of the sources which cover these topics will be 
necessarily abridged. This section is a survey of some, but by no means all, of the relevant 
literature regarding these topics. This material is mentioned in order to provide a context 
within which the methods employed are expected to operate. It is necessary to understand the 
limits of present analysis techniques in order to have realistic expectations for the use of those 
techniques in a given application. In this spirit, the following material was examined prior to 
and during the course of this investigation. 

Grid generation codes are employed to provide a convenient means of representing a 
complex geometrical shape in a way which can be modeled in a rectangular finite-difference 
grid. These codes typically transform the standard Cartesian coordinate system, which 
represents physical space, into a body-fitted curvilinear coordinate system, which represents the 
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computatuiona! space within which the governing equations are solved. Typically these 
methods result in the x-y coordinates corresponding to the grid points in the curvilinear 
coordinate system. The methods differ in the form of the transformation equation which must 
be solved and therefore in the boundary conditions which must be specified. Initial efforts at 
employing body-fitted coordinate systems were conducted by Winslow [12] in 1966 and Chu 
[13] in 1971. Both evaluated the Laplace equation as a means of transforming the coordinates 
from one system to the other. In 1973, Amsden and Hirt [l 4] used the same equations to 
provide a method of generating general curvilinear grids. Thompson, Thames, and Mastin 
[15-17] expanded this approach, during 1974-77, to include any number of bodies in the 
computational space. Their methods have been used quite extensively. In 1979, Steger and 
Sorenson [18] provided for angle and distance control of the grid lines at the inner boundary. 
In the following year, Sorenson [19] extended this method to control angles and spacing at the 
outer boundary. This method is employed in the so-called GRAPE code which is used for this 
investigation and will be described later. 

Other grid generation methods employ solutions of hyperbolic and parabolic partial 
differential equations and geometric techniques. Barth, Pulliam, and Buning [20] employed a 
hyperbolic grid generator, in 1985, to examine ‘exotic’ airfoils, including a glaze ice shape 
provided by the author. In 1978, Gibeling, Shamroth, aud Eiseman [21] developed a geometric 
grid generation technique, refined by Eiseman [22] which was used by Shamroth [23] in 1985 to 
evaluate steady and unsteady airfoil flow fields. 

Early efforts at evaluation of the Navier-Stokes equations considered incompressible 
laminar flow. Examples are those of Mehta and Lavan [24] and Lugt and Haussling [25], both 
from 1975. Mehta and Lavan solved these equations with a stream function-vorticity 
formulation to examine flow about an impulsively started airfoil. Lugt and Ilaussling also 
used a stream function-vorticity approach to examine flow about an abruptly started cylinder. 
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In an alternate approach, during 1977 Reddy and Thompson [26] applied an integro- 
differential formulation to the problem of incompressible flow in a doubly connected region. 
This was used to evaluate symmetric airfoils at zero angle of attack with a Reynolds number of 
less than 10 6 . Similarly, during the mid-1970’s, Wu and Sampath [27] and Wu, Sampath, and 
Sankar [28] applied an integro-differential formulation [29] to both an impulsively started 
airfoil and an oscillating airfoil. Finally, in 1980, Sugavanarm and Wu [30] attempted to use a 
two equation k-e turbulence model with a vorticity-velocity formulation. 

The primitive variable approach has also been employed for the evaluation of 
incompressible laminar flows. Harlow and Welch [3l] employed the Marker And Cell (MAC) 
method to investigate time dependent flow of a fluid with a free surface. This approach was 
developed further by Hirt and Harlow [32]. Later, Hodge [33] and Hodge and Stone [34] 
employed a successive over relaxation (SOR) iteration approach in a body fitted curvilinear 
coordinate system. The alternating direction implicit (ADI) approach was used by Ghia, 
Hankey, and Hodge [35] to calculate incompressible driven flow in a square cavity for Reynolds 
numbers under 1000. 

Compressible flow over airfoils has been and continues to be examined by a large 
number of investigators. Verhoff [36] applied MacCormack’s fully explicit method [37] to this 
problem but was restricted by small time steps in order to maintain numerical stability. 
Deiwert [38] also used this method to examine transonic flow. 

Implicit methods have been employed for the laminar compressible Navier-Stokes 
equations in an effort to avoid the stability limitations present in explicit schemes. Gibeling, 
Shamroth, and Eiseman [2l] applied the Briley-McDonald [39] formulation to examine 
dynamic stall. Sankar and Tassa [40] used a similar approach to examine an oscillating airfoil 
in a low Reynolds number flow. 


The examination of turbulent compressible flows presents an additional degree of 
complexity to the calculation of airfoil flow fields. The various investigators have generally 
employed either a zero-equation algebraic eddy viscosity model or a two-equation k-e model, 
where an n-equation model refers to the n-number of additional partial differential equations to 
be solved. Recently, there has also been considerable interest in the single ordinary differential 
equation model of Johnson and King [41]. All of these types of models have been compared by 
Coakley [42] for cases of shock induced separation, and the Johnson-King model was considered 
to have performed favorably. A k-e model that has been used by several invetigators (e.g, [23] 
and [43]) is described by Launder and Spalding [44]. The algebraic models most commonly 
used are those of Cebeci and Smith [45] and Baldwin and Lomax [46]. These models have been 
used sucessfully over a wide range of airfoil configurations for conditions below stall. At higher 
angles of attack, there have been some difficulties observed however. The models tend to 
overpredict the turbulence level and damp out some of the vorticity generation ocurring under 
these conditions. This is of particular concern for the iced airfoil condition due to the presence 
of the stationary separation bubble and its subsequent shedding at higher angles of attack. 
These considerations will be discussed in more detail in later chapters. 

Early investigation of the turbulent airfoil problem was performed in 1979 by 
Shamroth and Gibeling [47]. The method described therein employed a mixing length type 
turbulence model and was used to examine airfoils at low angles of attack. This approach was 
employed again during the following year by Shamroth and Gibeling [48] to examine airfoils in 
stall and again in 1981 by Shamroth [49] to examine airfoils pitching at low incidence. Tassa 
and Sankar [50] in 1981 and Sankar and Tang [51] in 1985 used an algebraic mixing length 
model to study dynamic stall. A k-e model was employed in 1985 by Shamroth [23] to 
investigate steady flow over a NACA0012 airfoil at low angles of attack. The basis for the 
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code employed in this investigation was developed by Steger [52] in 1978, who used the 
numerical scheme of Beam and Warming [53] and the Baldwin-Lomax turbulence model. 

The code developed at the Ames Research Center by Steger, ARC2D, has been 
further enhanced by Pulliam [54]. Pulliam introduced a diagonalization of the blocktridiagonal 
inversion for the implicit operators which resulted in a computationally more efficient set of 
scalar tridiagonal inversions along with a series of 4x4 multiplications. Additionally, he 
vectorized the code to take advantage of the capabilities of CRAY type computer architecture. 
This code has been used to evaluate a large number of airfoils over a large range of Reynolds 
number and Mach number conditions [55]. ARC2D and GRAPE will be the basic numerical 
tools used to examine the iced airfoil flowfield and will thus be described more fully in the next 
chapter. 
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CHAPTER 3 


DESCRIPTION OF COMPUTER CODES 

Solution of the Navier-Stokes equations with a finite difference method for a highly 
irregular geometry requires the use of two separate codes. One is used for the transformation 
of the physical space, which holds the object of interest, into a rectangular computational 
space, which is more suitable for the use of a finite difference code. The other is the actual 
Navier-Stokes solver. Both codes in this study were supplied by the NASA Ames Research 
Center. The grid transformation code is called GRAPE and was written by Sorenson [19]. 
GRAPE is an acronym for GRid transformation Algorithm for solution of the Poisson 
Equation. The Navier-Stokes solver is called ARC2D and was originally written by Steger [52] 
and later modified by Pulliam [54]. 

The following sections give a description of the basic equations being solved by the 
two codes used in this study. The form of the equations, as implemented in the codes, is also 
developed from these basic equations. Much of the development presented can be found in the 
references cited. Inclusion of this development is for the sake of familiarizing the reader with 
the specific mathematical formulation of the physical systems being modeled. This 
presentation is thorough enough to give the flavor of the calculation methods employed, while 
not purporting to be an exhaustive examination of work which has been described previously. 


13 



3.1 GRAPE Code Description 


The GRAPE code is based on a transformation of the physical x-y coordinate system to a 
body-fitted £-r) coordinate system through the use of a Poisson equation. The equations to be 
solved are, 


£xx + £yy — P 


(3.1a) 


Vxx + Vyy — Q 


(3.1b) 


where P and Q are constants which can be manipulated to control spacing and skewness of the 
resulting grid. Subscripts indicate differentiation with respect to the given coordinate. 

If the new coordinates are defined as functions of x and y then, 


t = «(*.y) 


(3.2a) 


and, 


v = v(x,y) 


(3.2b) 


It is desired to obtain the x and y coordinates of the rectangular £-17 computational 
grid. Thus, expressing (3.1) in terms of differentials of x and y with respect to £ and tj is 
necessary. The resulting partial differential equations are then solved using a computational 
technique described below. 

The relationship between (3.1) and differentials of x and y with respect to £ and rj is 
obtained by employing the Jacobian of the transformation from one coordinate system to the 
other, 


_ d(x,y) 

0(£,»?) 


(3.3) 


14 



or, 


J = X£ y f, - 

(3.4) 

Equation (3.4) can be rearranged to yield the following expressions, 


II 

X 

(3.5) 

r- 

X 

1 

II 

UJ* 

(3.6) 

vx — — 

(3.7) 

TJy = X^/J 

(3.8) 

Applying equations (3. 5-3. 8) to (3.1) yields, 


” 2/Jx^ + TXy W = “ j2 ( Px ^ + Q x *?) 

(3.9a) 

~ 2/?y^ + jy-qr] — — J 2 (Py^ + Qy^) 

(3.9b) 

where, 


2 , 2 
a — Xff + yrj 

(3.10) 

P = x^Xfj + y^yrj 

(3.11) 
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(3.12) 


In the GRAPE code, P, Q, and values of x and y are specified as input. The x and y 
values are input as a discrete set of points which can be used as is or redistributed using a 
curve-fit routine and a pre-defined distribution function. The boundary values correspond to 
the surface of the object being considered and an artificial boundary at some distance away 
from the object. The distance to the outer bondary is set at some value that is considered to 
be sufficient for free-stream conditions to prevail. The outer boundary can be specified as a C 
or 0 type grid. Equation (3.6) is solved by an iterative method employing a successive line 
over-relaxation solution procedure. The solution proceeds along lines that run in the £ 
direction. Convergence is obtained when the absolute value of the largest correction in x and y 
is below some desired level. Normally the convergence criteria is set as a drop in this 
correction value of six orders of magnitude. If convergence is not obtained, several relaxation 
parameters can be changed in an effort to improve the results. Also, the point distribution at 
the surface can be altered in an attempt to avoid regions of high curvature. The former 
approach is preferred, as this does not require alteration of the input geometry specification. 
Further explanation of the theory underlying the GRAPE code may be found in reference [19]. 

GRAPE has been used successfully to produce grids for many complex geometries, 
including ice shapes, as seen in figure 1. The output of this code is the x and y locations of the 
grid points from the uniform rectangular £-77 computational coordinate system. These 
locations are then used as input for the ARC2D code, specifying the locations at which the 
Navier-Stokes equations are to be solved. 

The ability of a particular grid to provide appropriate spatial resolution for a given 
flow-field is a question of much debate. Some of the influential factors are: the spacing of grid 
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Figure 1(a) Grid for rime ice leading edge 
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points near a solid surface in the normal direction, the spacing of grid points in the streamwise 
direction, the distance from the solid surface to the free-stream boundary, and the skewness of 
grid lines in regions of high curvature. Presently, the only measure of a grid quality is the 
degree to which the grid generation code has converged. Certainly, this allows for evaluation 
of how well the transformation equation has been modeled. However, the degree of 
convergence does not give an indication of whether the grid spacing is appropriate for 
resolution of the important flow phenomena. These concerns, important as they are for any 
computational fluid dynamics problem, are not the major emphasis of this research. As such, 
they will only be discussed in conjunction with specific problems that developed during the 
course of the research. 

The grids used in this work were similar, in degree of spatial resolution, to those 
employed for analysis of viscous flow over a clean airfoil (see, for example, Pulliam [54]). The 
spacing between the inner and outer boundaries was on the order of ten chord lengths. This is 
considered appropriate to provide free-stream conditions at the outer boundary [55]. The 
spacing of the first grid point off of the airfoil surface was on the order of .00002 of a chord 
length. This is sufficient spacing to provide approximately twenty grid points in the airfoil 
boundary layer. Grid spacing around the surface was concentrated near the leading edge in 
order to resolve the geometry of the ice shape. This is accomplished by providing a greater 
number of input points in this region and altering the point distribution function mentioned 
earlier. The high curvature of the ice shape results in significant skewing of the grid in this 
region. Concentration of points near the leading edge also results in scarcity of grid lines in the 
wake region. This can result in improper representation of the wake and must be considered 
carefully when creating a grid. 
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3.2 ARC2D Code Description 

The code selected for solving the Navier-Stokes equations, ARC2D, was originally 
developed at the NASA Ames Research Center, as mentioned previously. This code can be 
configured to solve the thin-layer Navier-Stokes equations or it can additionally include the 
explicit £ and cross derivative terms. 

The full equations are shown below in Cartesian coordinates. 


where, 


and, 


with, 


Q t + E x + Fy = Re'^Evx + Fv y ) 


(3.13) 


Q = 


p 


P u 


pv 

pu 

py 

, E = 

pu 2 + P 
/?uv 

, F = 

puv 

pv 2 + P 

e 


u(e + P) 


v(e + P) 


( 3 - 14 ) 



0 


0 

Ev = 

r xx 

II 

> 

r xy 


r xy 


r yy 


e 4 


U 


(3.15) 


r xx — j (4u x — 2 Vy ) 


(3.16a) 


19 



T xy = ^(“y + v x) 


(3.16b) 


Tyy = | (— 2u x + 4Vy) 

(3.16c) 

e 4 = ur xx + vr xy + H Pr'^T - 1) 1 ( ft2 )x 

(3.16d) 

f 4 = ur xy + vryy + n Pr _1 (T ~ 1 )’ 1 ( a2 )y 

(3.16e) 

P = (7 - l)^ e “ \ Pi" 2 + y2 )) 

(3.16f) 

These equations are then transformed to the body-fitted coordinate system 
established by the grid points taken from the GRAPE code output. According to the 
development presented by Vinokur [56], the strong conservation law form of (3.2) can be 


maintained for new independent variables of the form, 


€ = £(x,y,t) , t) = »?(x,y,t) , t = t 


Retaining strong conservation law form is important in that it is possible to 
difference the equations by a variety of stable schemes, each of which can be chosen so as to 
conserve mass, momentum, and energy for the total flow region. The time dependence is 
shown for these transformations and was not included in the GRAPE code formulation. Since 
only the x-y positions of the ^ grid are transferred from GRAPE to ARC2D, this dependence 
is not necessary for the GRAPE code. ARC2D uses the x-y coordinate information to form the 
metrics again internal to the code. The roles of the independent variables are reversed for the 
ARC2D code and hence the Jacobian has the form, 
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j ' 1 = *£ y Tf ~ x tj 


(3.17) 


The contravariant velocities are defined along the £ and 77 coordinates as, 


U — + £x u + £y v 

(3.18) 

V = T) t + t] x U + Tj y v 

(3.19) 


where the metrics are formed by 
the two coordinate systems. The 
coordinate systems. 


considering the differentials of the independent variables of 
metrics describe the following relationships between the two 


£x — Jy 77 1 £y — — J x 77 > “ ~ x r€x ~ YtZ y 


lx — Jy^ > Vy — “ > Vi — ~~ x tVx “ YrVy 


(3.20) 


Applying these transformations to the governing equations (3.2) results in the 
following form, 


Q r + ]L + F, = Re'^Ev^ + Fv,,) (3.21) 

where, 



P 


P u 


P V 

Q = J' 1 

pu 

pv 

' , E = J' 1 

puU + £ X P 

/>vU + £yP 

, F = J’ 1 

puV + 7?x p 

PVV + TJyP 


e 


u(e + P) - * t P 


v (e + P) - T) t P 
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Ev = J ^Ev^x + Fv£y) , Fv = J ^Evr/x + Fvrjy) 


(3.23) 


The thin-layer form of the equations is obtained by neglecting gradients of the 
viscous terms in the streamwise direction (i.e. the £ direction). This is similar to a boundary 
layer type assumption, however, unlike the boundary layer equations, the cross-stream pressure 
gradient is retained. This allows evaluation of regions with recirculation. 

Selection of full or thin-layer equations is dependent on the phenomena being 
modeled and the suitability of the grid. If the grid does not have sufficiently fine spacing in 
the streamwise direction, then use of the full Navier-Stokes equations may not be warranted. 
For airfoils with attached boundary layers, the increased run times required for the full 
equations are not justified by a significant increase in accuracy. This may not be the case for a 
separated flow field. Therefore, due to the leading edge separation which can occur with an 
iced airfoil geometry, the solution of the full Navier-Stokes equations should be examined. The 
thin layer equations have the form. 


Q t d* E^ + F Tj — Re Fv,j 


(3.24) 


Employing equations (3.15) and (3.23), the right hand side of (3.24) can be rewritten 


in the following form, 


Re _1 Fv^ = Re' 1 (j'^Evtjx + Fvijy)^ 
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or, 


0 


Re^F Vjj = Re' 1 J' 1 


T xxVx + r xy 7 7y 
T xy T lx + r yy 7 ?y 

e 4 T Jx + ^'Hy 


where Txx, r X y, r yy> e 4 > an< i ^4 are gi ven in (3.16). The shear stress terms can be 
with the u and v derivatives expanded by the chain rule, 


r xx ~ (A 4* 2/i)(£ x ii£ + ^x u t;) + ^(£y v £+ 7 7y v r;) 
r xy = /*((£yu^ + Vyu^) + (£ x v^ + Vx^rj)) 

Tyy = (A + 2/i)(£yV£ + TJyV^) + A(^ X U^+ Vx^l]) 

and the transformed e 4 and f 4 terms become, 


e 4 = U1_ xx + vr xy + ^Pr _1 (7 - ^(txd + T) x d, 1 a 2 ') 
U = u^xy + vr yy + /iPr _1 (7 - l)' 1 ^ y ^a 2 + 77y^a 2 j 


Substituting (3.27-29) into (3.26) yields the following relations, 
Re' 1 J" 1 ^Ev^ x + Fvqyjq = Re _1 S^ 


(3.26) 
rewritten 

(3.27) 

(3.28) 

(3.29) 

(3.30) 

(3.31) 

(3.32) 
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where, 


/j(t7x 2 + 7 7y 2 ) u V + (a*/3)jJx(»Jx u »/ + Vy^rj) 


S = J' 


^('/x 2 + »7y 2 ) u ?7 + (/ i /3)f?x(’7x u »j + *?y v r/) 


«Pr ' 1 (7 — l)' l (r/ x 2 + , ?y 2 )( a Tj) 2 +(/ J /2)( , 7x 2 + *7y )( u + v )r? 
+ (a</ 6 )( , ?x 2 ( u2 )j? + r ?y 2 ( v2 )»7 + 2»7x , /y( uv )j;) 


(3.33) 


Finally (3.24) and (3.25) are rewritten in a more compact form, 


Qr + + F jj — Re Sjj 


(3.34) 


This form of the Navier-Stokes equations can now be solved by forming an equivalent 
finite difference representation. Applying the implicit three point time differencing scheme of 
Beam and Warming [53] yields an expression of the form, 


aq" = IM (AQ")t + + J~T4 A 4"" 

+ 0 - \ - <t>) At 2 + At 3 


(3.35) 


where AQ n = Q n+1 — Q n and Q n = Q(nAt). The values of 0 and <j> are chosen to provide 
either first or second order accuracy in time. Typically, when the code is run in a time 
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accurate mode for unsteady calculations 6 is set to 1 and <f> is set to 0. This results in second 
order time accuracy for these calculations. For steady-state calculations, a spatially varying 
time step is employed to speed convergence. In these cases, the values of 0 and 4> arc set to 1 
and 0.5 respectively. This yields only first order accuracy in time. This is acceptable as long 
as there is convergence to a steady-state solution. 

The equations are further modified by employing a local time linearization as 
described by Pulliam [54], The equations then take the delta form of the algorithm, 


I + h ^A n ^ + h 


Re" 1 



= - h 


(e") { + (F")„ 



(3.36) 


where A = <9E/<9Q, B = dY/dQ, and M = 5S/5Q. The right hand side of (3.36) is the 
explicit part and the left hand side of (3.36) is the implicit part of the algorithm. 

The form presented in (3.36) is for the thin-layer Navier-Stokes equations. As noted 
before, ARC2D has the capability to include the streamwise and cross-term derivatives as an 
alternative to the thin-layer approximation. This is accomplished by retaining the explicit 
portions of Ev in the algorithm. The use of these terms was examined with respect to the iced 
airfoil flow-field and found to have no influence on results for the grid systems used. 

The spatial differencing employed by ARC2D is second order central differencing. 
Upwinding is also included as an optional approach but is generally used only for resolution of 
shocks. The matrix resulting from application of the central differencing is a 
(Jmax * Kmax * 4) x (Jmax * Kmax *4) square banded matrix. This form is sparse but very 
computationally expensive. In order to decrease the run times, the solution process is 
simplified by an approximate factorization [53] of the two dimensional operator (3.36) into two 
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one-dimensional operators of the form, 


[ I + h (A n ) J [ I + h (b"),, - h AQ n 


(3.37) 


= - h 


(E n )^ + (F n )»7 “ Re 


The solution algorithm now consists of two implicit operators each of which is block 
tridiagonal. The program flow now consists of two one-dimensional sweeps, one in the £- 
direction and one in the redirection. The resulting procedure is more economical both in terms 
of run time and computer storage [54]. 

Artificial dissipation is added to the algorithm to ensure stability, especially for 
transonic flow with shocks. This is required because, in high Reynolds number viscous flows, 
scales of motion exist which cannot be resolved by the numerical scheme. In the actual 
physical problem, these high frequency waves are brought about by the interaction of the 
convective terms in the momentum equations. This is normally accounted for by viscous 
dissipation but, since the scales are not resolved by the code there is no mechanism for 
removing this energy from the solution. 

ARC2D employs two techniques for introduction of artificial dissipation. Since 
ARC2D can be run as either a steady-state or unsteady code, two forms of time differencing 
are used as mentioned earlier. For the unsteady time-accurate mode, explicit fourth order and 
implicit second order smoothing with constant coefficients is included. Typically the fourth 
order explicit dissipation is set to zero for fine grids, such as those used in this investigation. 
In steady-state calculations, nonlinear artificial dissipation of mixed second and fourth order is 
employed. This dissipation model is based on the work of Jameson et al. [57]. A more 
detailed description of these models and their effect on convergence is provided by Pulliam [54]. 
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3.3 The Baldwin-Lomax Turbulence Model 


Theoretically, the Navier-Stokes equations fully describe the dynamics of a viscous 
flow, whether laminar or turbulent. In practice however, it is not feasible to evaluate the 
finite-difference equations on a grid fine enough to resolve all the details of the turbulent 
structures. This implies the need for some form of turbulence model. Typically, algebraic, and 
two-equation models are presently being employed in aerodynamic codes. ARC2D employs an 
algebraic model developed by Baldwin and Lomax [46], This model is a variation of the model 
developed by Cebeci and Smith [45] for boundary layer analysis. The Baldwin-Lomax model 
was developed to avoid the need for determination of the displacement thickness, which is a 
natural diagnostic of the Cebeci-Smith method but difficult to evaluate with a Navier-Stokes 
code such as ARC2D. This is due to the difficulty in defining the local free-stream velocity on 
a non-rectangular grid. 

The Baldwin-Lomax model divides the turbulent boundary layer into an inner and 
outer region. The eddy viscosity, \i u is then evaluated by examination of the vorticity level in 
each region. The inner region encompasses the laminar sublayer and the buffer region. The 
outer region is the wake-like region of the boundary layer and hence is modeled by a locally 
constant length scale and a velocity scale dependent on the vorticity. 

The form of the eddy viscosity equation in the inner region is, by dimensional 
reasoning, proportional to the density times a typical length scale times a typical velocity 
scale. 


f‘t = P * v inne 


The length scale is given by, 

£ = K y [l - exp(-y + /A + )] 


(3.38) 


(3.39) 
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where k is the von Karman constant, y is the normal distance from the wall, A is the van 
Driest damping constant and, 


+ w m 

y = y — jr- w 


(3.40) 


The velocity scale is taken to be v inner = £ \u\ ) hence 


n t = P i 2 M 


(3.41) 


In the outer region, the eddy viscosity is given by, 


Vt — K C C p p ¥ wake Fjciebiy) 


(3.42) 


where K is the so called Clauser constant, C C p is an additional empirical constant used by 
Baldwin-Lomax, F He6 (y) is a factor which tries to account for intermittancy and will be 
described later. F wake is a function defined by Baldwin and Lomax as, 


= Min 


y max Fmaj 


(CwJc U difj 2 ymaxj/Fmaa 


(3.43) 


where F ma x is the maximum value of the function, 


F(y) = y M [l - exp(-y + /A + )J 


(3.44) 


along a grid line in the ^-direction. The y-value at which F(y) reaches this maximum is 
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designated as y max- The value of U di j j is given by, 

U d , 7 / = (u 2 + max - (u 2 + V 2 ^ 2 mi „ (3.45) 

This formulation results in the length scale being determined by the distribution of vorticity 

along a normal to the surface. 

The inner region model extends outward until the value of fi t obtained with this 
model equals the / i t value obtained with the outer region model. The interface value is then 
taken as the average of the two model values at this point. The outer region extends from the 

point of equality with the inner region to a distance at which the velocity is equal to the local 

free stream value, which is taken to be the value at the outer boundary of the solution domain. 

The effect of intermittency has been modeled by incorporation of a variation of the 
Klebanoff intermittency factor. This term is given by, 

r 6 -|— 1 

F«,*(y) = 1+5.5 (3.4C) 

where 0.3. The value of F kUb is essentially unity for small values of y and drops to 

almost zero for large values. The transition from unity to zero occurs rapidly at y values close 
to y max- 

The Baldwin-Lomax model works very well for attached boundary layers and for 
small separated regions. Large separated regions can cause some difficulties for the model 
which results in improper representation of the turbulence in the region. In an attempt to 
rectify this problem, an alternate turbulence model has been developed which employs an 
approach similar to the Baldwin-Lomax model but which tries to avoid some of the difficulties 
encountered in the large separation regions. This is also a Modified Mixing Length (MML) 


29 


model, which bases the determination of the mixing length on the wall shear and the distance 
from the wall until some maximum is reached. The value of this maximum length scale is 
determined by evaluation of an attached boundary layer in a parametric study. The velocity 
scale is determined by the vorticity level and the mixing length at each point in the flow. This 
approach avoids the problem of trying to determine which of the local maxima of the Baldwin- 
Lomax F-function should be used to evaluate the viscous region on the iced airfoil. In the 
following chapters, the MML model will be described and compared to the Baldwin-Lomax 
model, especially in regard to calculation of cases at and above stall. 


i 
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CHAPTER 4 


MODIFICATIONS FOR ICING ANALYSIS 

This chapter details the changes made to the ARC2D code in order to enhance its 
ability to calculate the iced airfoil flowfield. Specifically, the approach to modeling 
turbulence was replaced by a simple but workable model in an attempt to provide a more 
realistic distribution of turbulent viscosity throughout the separated flowfield that develops 
behind the horns of the ice shape. The modifications were made in order to calculate the 
separated flow aft of the horns and to allow evaluation of the premature stall of the iced 
airfoil. 

Navier-Stokes codes have been used to evaluate separated flowfields (e.g. Mehta and 
Lavan [24]), but successful attempts have been confined to cases with moderate Reynolds 
numbers. For higher Reynolds numbers, typical of a turbulent flow, the ability to predict the 
maximum lift has not been demonstrated. This indicates that some attention to the 
turbulence model is necessary for determination of stall in an airfoil under normal operating 
conditions. The special circumstances of an airfoil with a leading edge ice accretion allows 
definite identification of the stall mechanism. An attached recirculation region aft of the 
horns breaks away from the surface which results in massive separation and stall. This 
eliminates one of the difficulties experienced in evaluation of stall on a clean airfoil. Thus, 
turbulence model alterations can be directed toward better evaluation of a specific flow rather 
than for employment in a more general and therefore more widely varying flowfield. 
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4.1 The Baldwin-Lomax model examined 


! 


The turbulence model originally employed in the ARC2D code is an algebraic eddy 
viscosity model developed by Baldwin and Lomax [46] and was described in some detail in 
the previous chapter. Calculations using this model for an attached boundary layer have 
been very successful [54]. However, the use of this model in separated boundary layers has 
previously resulted in some difficulties, as described by Degani and Schiff [58]. In their case, 
evaluation of cross-flow separation on a pointed cylinder at high angle of attack, the model 
suppressed secondary vortices when used as originally implemented. The difficulties 
experienced were attributed to the evaluation of the maximum of the F-function, described in 
Eq. (3.44). 

The presence of multiple maxima in the F-function led to excessively large eddy 
viscosity values. This in turn damped out the smaller flow structures, including the 
secondary vortices, and resulted in the center of the main vortex being displaced from the 
surface. Degani and Schiff modified the Baldwin-Lomax model by essentially deciding a 
priori which maxima was the more significant for their calculation. Their results indicated a 
significant improvement in determining the details of the flow T separation for the geometry 
being considered. 

Similar problems were experienced in this study during evaluation of the iced airfoil 
at high angles of attack. In this case, the F-function developed a third peak, as shown in 
figure 2. This type of profile was due to the large value of the vorticity at the wall and the 
smaller, although significant, vorticity levels in the vortex being shed from the surface of the 
airfoil. Due to the transient nature of the flow field, this type of profile moves along the 
surface with the vortex, further complicating the selection of the appropriate value of Fmar. 

As a result, the turbulence levels in the flow on the upper surface have steep 
gradients and change along with passage of the vortex. This is shown in figures 3 and 4, 
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which indicate the pressure distribution, stream function contours, and eddy viscosity levels 
at different points in time for an iced airfoil at an angle of attack above stall. As shown in 
the figures, the eddy viscosity levels increase and decrease abruptly in regions of high vorticity 
associated with the vortex travelling along the airfoil surface. Also apparent is the fact that 
the eddy viscosity levels are tied to the motion of this vortex and that points outside this 
region have little, if any, turbulent dissipation. It would seem more appropriate that some 
turbulence would remain in the regions between vortices and that the sharp gradients seen in 
the figures would not be present in an actual turbulent flow. 

The distribution and overall level of the eddy viscosity can result in significant 
changes in the size and shape of vortex patterns in the flowfield, as indicated by the 
experience of Degani and Schiff. The size and shape of the vortices and their shedding 
frequency all contribute to the time-averaged pressure distribution experienced by the stalled 
airfoil. Thus, it is reasonable to assume that by altering the development of this vortex 
shedding mechanism, the pressure distribution and hence the lift of the airfoil can be 
adjusted. If this adjustment is approached in a rational manner, then perhaps the lift of the 
stalled airfoil can be determined with a higher degree of accuracy than with the present 
turbulence model. 

With the complications of the Baldwin-Lomax model in mind and the goal of a 
more attractive model for separated flows, an alternate turbulence model was developed for 
evaluation of the iced airfoil. When selecting a new turbulence model it is appropriate to 
consider the degee of detail required of the model, the complexity of the model, the difficulty 
in implementing the model, and the resulting benefits of employing a given model. For 
aerodynamic codes the nature of ARC2D, the choice seems to lie between algebraic models 
and two equation models. Presently an effort is underway [59] to incorporate both a two 
equation model, that of Gorski [60], and the ordinary differential equation model of Johnson 


36 



and King [41] into the ARC2D code for evaluation of separated flows. Thus, in order to 
avoid an overlap of effort and to respond to the special circumstances of the iced airfoil 
conditions, an alternate approach was employed in this investigation. The model developed is 
essentially a mixing length model which seeks to avoid selection of the appropriate length 
scale based on the local maxima of an ad hoc function. 

4.2 The Modified Mixing Length Model 

The Modified Mixing Length (MML) model is a zero equation model based on the 
Prandtl mixing length hypothesis. That is to say that there are scales of motion for a 
turbulent flow associated with the transport of momentum much the same as the mean free 
path of kinetic theory characterizes the motion of gas molecules. These scales, if they can be 
determined, can be used to characterize the time averaged or mean flow motion of the fluid. 
The difficulty, of course, lies in the determination of the appropriate scales. 

The MML model was developed in an attempt to model the appropriate scales of 
motion for a flowfield which is alternatively separated from and attached to an arbitrary 
geometric surface. Observation of the motion produced by this surface, as shown in figures 3 
and 4, leads to some assumptions as to how the turbulence may be generated and transported 
in the flowfield. These assessments of the nature of the turbulent flowfield are then 
translated into a numerical mechanism for the introduction of turbulent viscosity into the 
calculation. This process requires some prior evaluation of the flowfield of interest as well as 
some reference to the work of other investigators which may be relevant. As such, some 
background for the development of this model will be presented next. 

The MML model tries to address the problem of distribution of turbulent 
dissipation throughout a separated flowfield. In so doing, it is necessary to have some 
indication of the character of the flowfield being examined. An iced airfoil is typified by 
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regions of flow separation aft of the protruding horns (see Figure lb)* These regions increase 
and decrease in size with varying angle of attack (AOA). The upper surface separation region 
increases with increasing AOA, while the lower surface region decreases in size. The opposite 
occurs with decreasing AOA. At some critical AOA, the bubble no longer reattaches to the 
surface and a large unsteady vortex shedding pattern develops. These two different flow 
patterns are shown in figures 5a and 5b, respectively. 

The closed separation bubble of figure 5a is modeled adequately by the Baldwin- 
Lomax model, as will be discussed in the next chapter. The unsteady flow pattern of figure 5b 
: is the flow which resulted in the difficulties discussed in section 4.1. Therefore, the special 

circumstances of this flow were examined and modeled during the development of the MML 
! model. 

5 The MML model is based on the expression associated with the Prandtl mixing- 

; length theory [61]. Basically the turbulent viscosity is taken to be, 

; Pt = PP M ( 4 -!) 

i 

i 

I 

The question is what to use for the evaluation of the mixing-length, L The mixing- 
length is dependent on distance from the wall. In an attached boundary layer theie are 
typically three regions; 

t ~ y 2 y + < 10 (4.2) 

I , 

t ~ ny 10 < y < 100 (4.3) 

m 

£ = const. y + > 100 (4-4) 

s 

i 

9 

I 
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where, 


y 


+ _ yjLr 

~ V 


( 4 . 5 ) 


and 


1 



( 4 . 6 ) 


and 


r w = 



( 4 . 7 ) 


The flowfields which seem to develop on the iced airfoil above stall have three 
distinct regions. Very near the wall is a region of low velocity and high vorticity, 
distinguished from conventional attached boundary layers only by the far field boundary 
condition (i.e. the vortex structures in the separated region). Further out is the moderate 
velocity, moderate vorticity recirculation regions of the main separated flow. Above this is 
the high velocity, low vorticity region of the outer flow. It would seem then that the 
vorticity is generated at the horns and along the surface of the airfoil. This observation is in 
line with the expression for development of vorticity in a Newtonian fluid given by, 


^ = (fi-V)U + 



Vx 


V x 


1 

P 


V x 


(pQ) 


+ v 


( 4 . 8 ) 


which is obtained by taking the curl of both sides of the equation of motion in vector form. 
In this expression, Q is the vorticity vector, U is the velocity vector, P is the pressure, p is 
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the fluid density, /z is the absolute viscosity, and ^ is the substantial derivative. The fust 
term in this equation indicates the change in vorticity in the flowfield due to strain rate. The 
second term is a source of vorticity due to the pressure force. The final term is the diffusion 
of vorticity by viscosity. In a two-dimensional flow the first term disappears and only the 
source and diffusion terms remain. 

As a solid body passes through a viscous fluid with initially no vorticity, the 
pressure gradient at the surface produces vorticity tangential to the surface. This vorticity is 
then spread through the flow due to the viscous dissipation. The linkage between this 
generation and diffusion of vorticity and the level of turbulence is provided by (4.1). 1 hus, if 
the vorticity is generated near the surface, it may be reasonable to assume that the value of 
fi t is associated with the flow near the surface. Also, as the vorticity is spread through the 
flowfield, the turbulence level develops along with it, hence the dependence of on |w|. 

Equation (4.8) indicates that there is no vorticity production in the separated region 
due to the absence of a significant source term. Therefore, it is reasonable to conclude that 
there is no further enhancement of the turbulent viscosity. This is embodied in Prandtl’s 
observation that the length scale very far from the surface is a constant, Eq. (4.4). This 
allows the turbulent viscosity to diminish along with the vorticity as would be expected from 
equation (4.8). 

The MML model is based on the idea that the length scale is dependent on 
conditions near the surface and that its level remains constant in the separated region. The 
length scales are thus established by conditions at the surface and are then transported into 
the separated flow regions along with the mean flow. The ultimate level of the turbulent 
viscosity in these regions is established by the length scale and the level of vorticity. This 
leads to a two layer type model as employed by Cebeci-Smith and Baldwin-Lomax, with the 
cross over point being established by conditions near the surface. 


41 


The inner region of the MML model encompasses the laminar sublayer and the 
logarithmic buffer layer. There are several empirical formulas available which are used to 
evaluate this region. In this case, the van Driest formulation is used and is given by, 


f(y) = «y 


(-y + /A + )' 


(4.9) 


The outer region is based on the observation that the length scale in an attached 
boundary layer saturates at a level of about 0.10(5, as shown in figure 6 taken from reference 
[62]. For a separated flow, there is no definite boundary layer thickness, hence the length 
scale in this region is defined with respect to the value of y* where, 


y 


* 




w 


The length scale is simply, 


(4.10) 


£ = const, x y* 


(4.11) 


with the constant to be defined empirically. 

The two regions are blended into each other through a function of the form, 



(4.12) 


where C^y* is the distance above the surface at which l saturates and C 2 controls the 
curvature of the blending region. The form of the mixing length profile as a function of 
distance from the surface is shown in figure 7. Note the similarity to the i curve in figure 6. 
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The MML model is thus defined by the following relations, 


y + < Cj 


^(y) = K y* 




(-y + / A+ )' 


(4.13) 


and. 


y + > C x <(y) = « gi y* 


where 


+ 


can be rewritten as, 


(4.14) 


y 



(4.15) 


From (4.14) we see that the constant in (4.11) is given by, 


* C i 
const. = K pr 


(4.16) 


where the values of C l and C 2 can be varied to match empirical results. 

This model is dependent on the value of r w and thus will in general enhance the 
level of turbulent viscosity near regions of separation and reattachment. On the other hand, 
the value of r w in the backflow region of a separated flow is relatively large, resulting in 
lower values of the turbulent viscosity. This agrees with a number of the observations of 
Simpson et.al. [63] regarding two-dimensional separated flows. These are that the separating 
shear layer behaves progressively more like a free shear mixing layer and that the part of the 
backflow adjacent to the surface has little Reynolds shear stress effects. There are some 
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observations which are not incorporated however, but these appear to be more characteristic 
of a closed separation bubble than of the massive separation being considered here. 

As a result of the use of r w near regions of separation and reattachment, there is a 
possibility of the mixing length becoming excessively large. This is seen in figure 8 a which 

shows |r w | 2 , hence the values of y* and £, near such a region. This problem is avoided by 

the use of a local average for r w - In the model used for this investigation, the r w is filtered 
spatially using the expression, 


r i,l “ °' lr i-2,l + 0-2r i-l,l + °' 4r i,l 
+ °- 2r i+l,l + 01r i+2,l 


(4.17) 


where the subscripts indicate grid points in the £ and 77 directions respectively. The 77 - 
direction subscript is set at one to indicate values at the surface. 

This spatial filtering assures the use of a non-zero value in the denominator of y* by 
converting a profile similar to figure 8 a to one more like 8 b. The use of this averaging also 
reflects the fact that separation and reattachment are processes extending over a certain 
region rather than an isolated event. Hence, it is expected that the MML model will capture 
the physical characteristics of these processes. 


4.3 Evaluation of the MML model 

The MML model was developed to address some of the discrepencies resulting from 
use of the Baldwin-Lomax model with separated flows. The new model should also perform 
reasonably for attached flows, since it incorporates all the elements of mixing-length models 
used for attached flows. In order to evaluate this model, a number of comparisons were made 
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between results from ARC2D using the Baldwin-Lomax model and the MML model for 
attached and separated flows on a NACA0012 airfoil with no ice shape. Comparisons to 
experimental results are also made when the latter are available. 

The cases examined are 0° AOA and 15° AOA for a Mach number of 0.12 and a 
Reynolds number of 1*41x1 0 6 . The grid used is shown in figure 9. There are 253 nodes in 
the ^-direction and 64 nodes in the 77-direction. There are 42 points along the wake cut and 
thus 211 points along the surface of the airfoil. The spacing of the first grid point normal to 
the airfoil surface is 2xl0 " 5 chord lengths. This grid is similar in size and spacing to the 
grids used for the iced airfoil. 

0° AOA-Attached flow 

The attached flow at this angle of incidence is a steady flow condition and thus the 
code can use the form of the algorithm designed for such flows. The time step used is thus a 
spatially varying value which results in the optimum convergence at each point in the grid. 
The value of the time step input to the code is thus a parameter used in an expression such 
as, 


At o 



(4.18) 


where At 0 is the input time step and is normally chosen to be 0(1). For the cases run here, 
the value of At 0 was set at 0.9. 

In evaluating the MML model, the values of C 2 and C 2 can be manipulated in order 
to produce results which agree with experimental information. By employing this process for 
the clean airfoil it is expected that the model will then be ’tuned’ and ready for use with the 
iced airfoil. The MML model was thus used with several values of Cj and C 2 and compared 
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(b) 

Grid system for NACA0012 airfoil a) Complete grid (253x64) b) Close-up of grid 
near airfoil surface. 


9 




to the Baldwin-Lomax results and to experimental information. The values used for the 
MML model are selected such that they produce length scales similar to those near the region 
of transition from inner layer to outer layer as shown in Bradshaw [64]. Variation of the C 2 
parameter produced overall effects similar to those of varying the C x parameter and results 
can therefore be reported for the C x parameter alone. The results of this comparison for drag 
values are shown in Table L 

These results indicate that the MML model can have significant variations in drag 
as a result of varying C x . However, the Baldwin-Lomax model produces drag values closer to 
experimental results. This discrepancy is examined further by looking at velocity gradients at 
the surface of the airfoil. The results for the Baldwin-Lomax model and the MML model, 
with the same C x values used previously, are compared in Table 2. 


Evaluation Method 

Cx 

c 2 

C D 

Baldwin-Lomax 

NA 

NA 

0.014 

MML 

300 

5.0 

0.019 

MML 

1000 

5.0 

0.020 

MML 

8000 

5.0 

0.016 

Experiment 

NA 

NA 

0.009 


1 

| 

Table 1 Force coefficients for several values of turbulence model parameters 
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Turbulence Model 

C, 

X 

C 

<9.1 

dy 

Baldwin-Lomax 

NA 

0.6 

2.0Sxl0 r> 

MML 

300 

0.6 

3.67 x 10 5 

MML 

1000 

0.6 

3.67x10 s 

MML 

8000 

0.6 

3.78 xlO 5 


Table 2 Velocity gradients at the surface as calculated using the Baldwin- Lomax 

and MML turbulence models 

This table indicates that the velocity gradients for the MML model arc 
approximately twice as large as the values from the Baldwin-Lomax model. If the total 
velocity profile at a given location is examined, it is found that the differences are relatively 
minor, the relative error based on u w being 0.01. Thus it can be seen that the values of the 
velocities very near the surface must be determined accurately to produce the correct drag. 
This is seen directly by examining the equation used to determine the frictional component of 
the drag. That is, 


Cj^ = C f (Ay sin a -f Ax cos a) 


where, 


ft; 


( d u _ dv\ 
y<9y dx j 

(booUoo 2 ) 


(4.19) 


(4.20) 
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The value of /i at the wall is essentially the molecular viscosity. The geometric 
values are the same for both calculations as are the free stream fluid properties. The only 
difference between the two calculations is due to the velocity gradients. Since the wall 

velocity gradients are calculated by using the first two grid points off the surface, the 

evaluation of these points is critical. These points correspond to a height of 7x10 5 chord 
lengths above the surface. Determining how the overall turbulence level affects the velocity 
gradient at the wall will require further study. Even in the case of the Baldwin-Lomax model 
the results differ from experiment by 52 percent. Abandonment of the MML model is not 
warranted strictly on the basis of these results. Indeed, as later results reveal, a strong case 
can be made for the adoption of the MML model. 

15° AO A - Separated flow 

The MML model was developed with this case in mind. At 15" AOA, the 
NACA0012 airfoil is near C T and under certain conditions, stall will occur. The 

L max 

NACA0012 airfoil has two possible modes for stall at this angle of attack. Gregory and 

O’Reilly [65] state that this airfoil stalls either from collapse of a short leading edge laminar 

separation bubble or from the rapid advance toward the leading edge of a trailing edge 
separation. The question of the conditions under which either mechanism occurs is not 
resolved by their experiment nor was it addressed in the work of Bragg [66]. 

As far as the calculations with ARC2D are concerned, there is no evidence of a short 
leading edge bubble using either turbulence model. There are differences, however, in the two 
calculations. The MML model produces a relatively large separation region near the trailing 
edge, while the Baldwin-Lomax model shows a much smaller separation region. This is shown 
in figures 10a and 10b, which show stream function contours for the Baldwin-Lomax model 
and MML model respectively. The difference between these two calculations is due to the 
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(b) 

Figure 10 Stream function contours for NACA0012 airfoil. AOA = 15\ 
a) Baldwin- Lomax model, b) MML model 
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distribution of turbulent viscosity near the trailing edge, as seen in figures 11a and lib. The 
Baldwin-Lomax model develops a large region of high fi t values which tend to suppress the 
development of the trailing edge separation. On the other hand, the MML model has high 
values of fi t only near the separation point, allowing the reverse flow region to develop 
downstream of this location. 

Figure 12 shows the pressure coefficients from these two calculations compared to 
the experimental values of Bragg [66]. As seen, the MML model captures the alteration of 
the pressure development at the trailing edge much better than the Baldwin-Lomax model. 
This is reflected in the values obtained with the two models. The Baldwin-Lomax model 
gives a C L value of 1.4 while the MML model gives a C L of 1.2. The experimental value is 
1.2. Thus, while the issue of the actual mechanism for stall is not settled, it can be seen that 
the MML model produces a more reasonable representation of the airfoil flowfield than the 
Baldwin-Lomax model for high angles of attack. 

Further testing 

Since the MML model was able to produce results which more closely model the 
actual flowfield for the stall condition, it was necessary to determine the appropriate values of 
C x and C 2 to employ. A series of cases near stall were run with different C x and C 2 values. 
The cases evaluated are indicated in Table 3. The values selected as the most appropriate 
from this examination will then be compared to similar results from an examination of the 
iced airfoil at stall. 
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Figure 12 Pressure coefFicient profiles for 15’ AOA. ARC2D with either the B-L model 

or the MML model versus experimental results of Bragg. 



Run Number 

AOA 

Ci 

^2 

1 

14 

1000 

5 

2 

15 

1000 

5 

3 

16 

1000 

5 

4 

14 

2000 

5 

5 

15 

2000 

5 

6 

16 

2000 

5 


Table 3 Evaluation of variation in on MML model results at several AOA 

values near stall 


These cases were all run in the time accurate mode. Thus, the lift histories can be 
examined to determine if there is some unsteadiness in the flow. The lift histories for runs 1- 
3 are shown in figures 13-15 respectively. These results indicate unsteady behavior starting at 
15° AOA. The time-averaged lift value is lower than the steady value at 14° and thus stall of 
the airfoil is indicated. The lift histories for runs 4-6 are shown in figures 16-18 respectively. 
These results indicate steady flow behavior and thus no stall of the airfoil. The resulting lift 
values indicate a progressive increase with AOA. 

The C l vs. AOA curve for a NACA0012 airfoil is shown in figure 19, taken from 
NACA TR-446 [67]. Stall of the airfoil is indicated at an AOA just above 16°. The 
value is 1.52 and the value drops off to 1.16 at 18". These results indicate that runs 4-6 
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Figure 16 Lift history for run 4 : AOA 14 , Cj — 2000, C 2 — 5. 
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Figure 17 Lift history for run 5 : AOA 15*, C 2 = 2000, C 2 = 5. 
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are a more accurate representation of the flowfield behavior. The value of Cj will be set to 
2000 due to better correlation with the experimental values. The C 2 value can be altered in a 
similar manner, with the additional effect that the velocity profiles will be altered. At this 
point, it is not an objective to match the velocity profiles exactly and the C 2 value will be left 
at 5. The velocity profiles will be discussed later along with considerations of transition 
location and grid refinement. 

The ARC2D code was next used to evaluate the NACA0012 airfoil at 18 AO A. 
This should be an unsteady result with significant vortex shedding. This unsteady behavior is 
manifested in the lower value shown in figure 19. The code was run employing the 
Baldwin-Lomax model, the MML model, and in an all laminar mode. The MML model used 
values of C x = 2000 and C 2 = 5, as prescribed above. The results examined are the lift 
histories, stream function contours, and eddy viscosity contours. Additionally the C j 
Op values are compared to the experimental information from NACA TR 446 [67} in Table 
4. The Cj and Cq values are found by time averaging over several periods of the shedding 
process. 


Evaluation Method 

C L 

C D 

S, 

Baldwin-Lomax 

1.65 

0.056 

NA 

MML 

1.35 

0.093 

0.03 

Laminar Flow 

1.08 

0.223 

0.16 

Experiment 

1.16 

0.19 

NA 


Table 4 Evaluation of turbulence modeling on post-stall pcrformacc prediction 
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Table 4 indicates that the MML model is better able to determine the Cj^ values for 
this post-stall behavior than the Baldwin-Lomax model. The MML model and the all 
laminar case results show reasonable differences from the experimental value, with the 
laminar case underpredicting and the MML model overpredicting the experimental value. 
The drag values show similar results. The Baldwin-Lomax model underpredicts drag 
significantly due to the attached flow behavior predicted using this model. The MML model 
underpredicts drag by an amount similar to the overprediction resulting from the laminar 
flow calculation. The laminar flow calculation results in a larger drag value due to the 
greater frequency of vortex shedding. It seems clear that the MML model performs better 
than the Baldwin-Lomax model for this case, but perhaps not any better than the use of no 
turbulence model altogether. The determining factor may be the evaluation of the unsteady 
vortex shedding phenomena. 

The lift histories are shown in figures 20-22. The vortex shedding is indicated by the 
periodic behavior of these curves. The shedding process can be characterized by the Strouhal 
number, S t , which is defined as, 


S, = (4.21) 

u oo 

where f is the shedding frequency, L is the characterisitic length, a is the airfoil angle of 
attack, and Uoo 18 the freestream velocity. For this case the characteristic length is the 
airfoil chord. The Strouhal numbers obtained from these lift histories are also shown in 
Table 4. 

The Strouhal numbers shown correspond to behavior that has been observed 
experimentally. The 0.16 value of the laminar flow case is approximately the value seen 
during bluff-body shedding. The 0.03 value of the MML model case has been documented in 
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Figure 22 


Lift history for NACA0012 airfoil calculation : AOA 18 , MML model. 



the work of Zaman and McKinzie [68]. The result obtained from the Baldwin-Lomax model 
seems to miss the unsteady stall behavior at this AOA entirely. The question of what is 
appropriate for use in further icing calculations is unclear at this point. A study of this 
behavior for an iced airfoil geometry is currently underway [69]. Results of that study could 
shed some light on whether the MML model is operating properly. 

Concluding Remarks 

The results of the examination of the clean NACA0012 airfoil with the two 
turbulence models suggests that each model has desirable characteristics. The 

Baldwin-Lomax model apparently produces drag values somewhat closer to experimental 
results for low values of AOA. The MML model however, allows modeling of the separated 
flow characteristic of a stalled airfoil. Certainly, it seems that further comparisons of the two 
models will be necessary. Evaluation of alternate well-documented airfoil shapes should help 
to indicate further the strengths and weaknesses of the two models. Particularly, the 
interaction of turbulence model and velocity gradient at the wall should be examined further. 
This would require extensive comparisons of the two models to measured values of velocity 
gradients, pressure coefficient distributions, and integrated force coefficents. 

This, however, is not the subject of the present study. For purposes of evaluation of 
iced airfoil performance, the pressing need is to determine the degradation of airfoil 
performance and to identify premature stall for a given ice shape profile. These needs 
translate into two requirements for the computational results. The code must be able to 
capture the main characteristics of the pre and post stall behavior for the iced airfoil. The 
pre-stall behavior includes determination of the size and shape of the attached separation 
bubble aft of the horns and the velocity gradients in that bubble. Additionally, the code 
must be able to determine the increase in drag and decrease in lift which can result from the 
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presence of the ice shape and the attached separation bubble. The code must also be able to 
determine the angle of attack at which the separation bubble detaches from the airfoil and 
the unsteady shedding process commences. The calculation of the value of also 

requires the code to predict the correct size, strength, and shedding frequency of the vortex 
which develops at the ice shape horn. 

Thus, the choice of turbulence model is still not clear at this point and the need for 
further testing is indicated. The ability of the MML model to capture the proper unsteady 
behavior near stall suggests that it may be possible to predict Cj using this model. The 
MML mode] also does not require selection of a maximum or minimum value of some 
variable from the mean flow field. This prevents the ambiguities found in the Baldwin- 
Lomax model for separated flow. The independence of this model from the solution 

procedure or grid structure should also allow easy transport of this model to alternate flow 
codes, including unstructured grid methods. The use of these two models will be further 
examined by using several iced airfoil geometries and will be discussed in the next chapter. 
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CHAPTER 5 


ICED AIRFOIL PERFORMANCE EVALUATION 

The tools described and developed in the previous chapters have been used to evalute 
the performance characteristics of several airfoil and ice shape combinations. The airfoils 
examined were chosen due to the availability of experimental data for verification purposes. 
The results for lift and drag coefficients using these airfoils have been calculated using both the 
Baldwin- Lomax model and MML model. Results indicate that a Navier-Stokes code coupled 
with the MML model can be used effectively to evaluate airfoil performance degradation due to 
icing. Calculations at angles of attack below stall indicate steady flow behavior with a 
recirculation bubble aft of the horns, when these are present in an ice accretion. At higher 
AOA values, the bubble separates completely and an unsteady vortex shedding process occurs. 
The MML model allows this process to develop and hence enables calculation of post-stall 
behavior. Comparisons between the Baldwin-Lomax and MML model for both pre and post 
stall behavior will be presented. The unsteady post-stall behavior, predicted using the MML 
model, will be examined in some detail. 

The needs of the icing community in regard to an evaluation of iced airfoil flowfields 
are twofold. Performance information is needed to evaluate the effects of ice accretion. This 
requires global integrated results such as lift, drag, and moment coefficients. Additionally, the 
particle trajectory codes and ice accretion codes require information on local velocities and 
pressures in the iced airfoil flowfields. These local results are also used to gain a greater 
understanding of the characteristics of iced airfoil aerodynamics. In that regard, stream 
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function contours, equi-vorticity contours, and eddy viscosity contours arc also useful for 
understanding the aerodynamic phenomena. The results have therefore been organized into the 
above two categories. For each airfoil/ice shape studied, both the global and local results will 
be examined in an effort to further understand the modeling capability of the code. 

The global results are the most prominent indication of the performance loss due to 
icing. The experimental results of Bragg, shown in figures 23 and 2 4, indicate the dramatic 
loss in lift and increase in drag produced by a glaze ice shape. Examination of these figures 
indicates that the lift diverges from the clean airfoil values only as the AO A approaches stall. 
The drag coefficient on the other hand is affected at all values of AO A. This points out the 
dual effect that ice accretions have on airfoil performance. Both the pressure induced forces 
and the frictional forces are altered by the presence of the ice accretions. The relative 
importance of each effect is dependent on the ice shape and the AO A. At low angles of attack, 
the pressure distribution is affected due to the separation bubbles, however the resulting 
changes to lift and moment are not significant. The changes to drag values are significant, and 
can be attributed to both pressure and frictional forces. At high angles of attack, the 
alteration to the pressure distribution is significant and leads to large changes in both 
quantities. 

These integrated force and moment coefficients are obtained by evaluating pressures 
and friction forces at each grid point on the surface of the airfoil. Both of these surface force 
values are then resolved into normal and tangential forces which in turn are transformed into 
lift and drag coefficients based on AO A. Monitoring of the development, with respect to 
iteration number, of these coefficients is one indication of convergence of the solution. These 
time histories of the force coefficients are also used to indicate whether a given solution is 
steady or unsteady. In the case of unsteady behavior, time-averaging of the pressure 
distribution is used to calculate the resulting forces on the airfoil surface. 
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Figure 23 Lift coefficient vs. AOA for a NACA0012 airfoil with artificial glaze ice shape. 
Experimental results, taken from Bragg [66]. 
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Figure 24 Drag coefficient vs. AOA for a NACA0012 airfoil with artificial glaze ice shape. 
Experimental results, taken from Bragg [66]. 
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The flowfleld associated with large glaze ice conditions, such as those shown in figures 
26 and 27, is characterized by a sizable disturbance to the flow just aft of the horns. This 
disturbance may or may not cause a significant alteration to the performance of the airfoil 
depending on the size of the accretion and the AOA of the airfoil. This is demonstrated most 
noticeably by the larger drop in Cj^ max reported by Bragg for the glaze ice profiles. Below the 
threshold at which stall begins, the alteration of the flow is restricted to the region 
immediately surrounding the ice accretion. This results in only minor changes to the lift as 
compared to the clean airfoil. The correct representation of the flowfield at these AOA s is 
important however from the standpoint of ice shape prediction and accurate determination of 
the drag. 

In order to provide code validation information for these lower AOA conditions, 
velocity and pressure measurements were taken by Bragg [66] for the NACA0012 airfoil and a 
glaze ice shape described below. The pressure measurements consisted of a densely packed 
series of pressure taps along the chord of the model, including the ice shape region and the 
entire upper and lower surfaces. This allowed very accurate evaluation of the C p distribution 
on the airfoil and of the pressure component of the force coefficients. The velocity 
measurements were taken with a split-film probe which was oriented in such a way as to allow 
determination of the magnitude and direction of the x-component of velocity. This was 
essential in evaluating the characteristics of the recirculating flow within the separation bubble 
aft of the glaze ice horns. 

Discussion of the flowfield characteristics at large values of AOA requires the 
definition of several terms which will be used frequently. The terms used will correspond to 
the definitions used by Mehta [24]. An attached separation bubble is a region bounded by a 
stream function contour of zero and the surface. The bifurcation point of the zero stream 
function contour is the separation point. The point of unification with the surface is the 
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reattachment point. A separation bubble is open or burst if it is not completely enclosed by a 
zero stream function contour and there is a closed contour within the region. A vortex is a 
region which is enclosed within equi-vorticity lines. Vortex shedding is the process of 
detachment of a vortex from the surface and subsequent convection into the free stream. 

The unsteady flowfield behavior at large AOA’s is the cause of the stall behavior 

reflected in the C T value and in the divergence of the curve. The three modes of stall; 

leading edge stall, trailing edge stall, and mixed leading-trailing edge stall; are described by 
Chang [69]. In the results described below, all three forms of stall have been observed. The 
ability of the code to indicate the type of stall and its ability to calculate correct force 
coefficients for these conditions, is evaluated thoroughly for the airfoil/icc shape combinations 
being considered in this study. A rime ice shape and two glaze ice shapes with two airfoil 
geometries are evaluated using ARC2D and the two turbulence models. 

5.1 Ice shape geometries 

The methods developed to evaluate iced airfoil performance must be independent of 

both the airfoil and ice shape geometry. This suggests that more than just one combination of 

airfoil and ice shape should be examined in order to increase confidence in the computational 
results. Unfortunately, there is not a large database of iced airfoil aerodynamic measurements 
available. The most complete information to date is that of Bragg [66] for the NACA0012 
airfoil. Bragg, Zaguli, and Gregorek [70] also measured the performance characteristics of a 
NACA63^-415 airfoil with rime and glaze ice shapes. The data from these two studies will be 
used to make comparisons with the computational results. 

The data available for the NACA63^-415 airfoil consists of pressure coefficient 
distributions at several angles of attack along with lift, drag, and moment coefficients for these 
same AOA’s. There was no investigation of the velocity profiles along the surface or in the 
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wake. The data available for the NACA0012 airfoil consists of both types of information just 
mentioned in addition to normal Reynolds stress ^ u /2 j profiles at several locations. The force 
coefficient values for the NACA0012 airfoil are available for AOA’s well past stall, while those 
for the NACA63^-415 are not. Thus, the global results will be presented for all three 
airfoil/ice-shape geometries while the local results will be concentrated on comparisons using 
the NACA0012 airfoil. 

Rime ice shape for the NACA63^-415 

The rime ice shape selected for use with the code corresponds to the R7 ice shape, 
described by Bragg et al. [70], and is shown in figure 25. This shape was selected for analysis 
because it appeared to have the profile least altered by the ice accretion. Indeed, at high 
AOA’s, this shape increased the lift of the airfoil over the range measured. This shape actually 
seemed to act as a leading edge flap. It was thus felt that this shape would provide an 
interesting test for the computational method. 

Glaze ice shape for the NACA63^-415 

The glaze ice shape selected for use with the code corresponds to the G3 ice shape, 
described by Bragg et al. [70], and is shown in figure 26. This shape was selected as a test for 
both the GRAPE code and the ARC2D code. The large concave region provided a critical test 
for the grid generation code. If this shape can be modeled accurately, then it is expected that 
the GRAPE code will be sufficient for the evaluation of most ice shapes. The large stagnation 
region between the horns and the rapid acceleration around these structures should also provide 
an important test of the ARC2D code. 
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Figure 26 Glaze ice shape for the NACA63^- 415 airfoil, (a) Close-up of leading edge, 
(b) Overall profile. 
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Glaze ice shape for the NACA0012 


This shape is the only ice geometry employed in Bragg’s later study [66] and 
represents a 5-minute glaze ice accretion. In this case, the geometry was a simulated ice shape 
and is shown in figure 27. This shape will be referred to as the G1 ice shape. The well-defined 
geometry, consisting of circular arcs and line segments, facillitates modeling in performance 
codes and allows for incorporation of the geometric detail deemed appropriate by the analyst. 
This geometry is the result of a coordination of effort between experimentalist and analyst 
typical of the icing program at NASA. 

5.2 Evaluation of Iced Airfoil Calculations 

This section examines the results of the ARC2D calculations for the ice shape 
geometries just defined. Comparisons are made to experimental information where available. 
The two turbulence models are compared in order to determine if there is any advantage 
obtained by use of the MML model for separated flow calculations. 

NACA63^-415, R7 ice shape 

This shape was tested in the NASA IRT wind tunnel and was run at a nominal 
Reynolds number of 5xl0 6 and a nominal Mach number of 0.14. The data were taken up to 
an AOA of 14.6 degrees. Comparisons between data and computations for Cj and C ^ arc 
shown in figures 28 and 29 respectively. Computed results are shown for both the Baldwin- 
Lomax and MML turbulence models. 

The lift calculations agree remarkably well with the experimental results up to an 
AOA of approximately 12.5°. The value of C^ was not determined experimentally for this 
shape. As seen in figure 28, the measured lift continues to increase over all AOA’s evaluated, 
with the exception of the value at 13.6" AOA. Examination of the pressure coefficient 
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Figure 29 Drag coefficient vs. AO A for NACA63^- 415 airfoil with R7 ice shape. 


Comparison of ARC2D results using both turbulence models to 
experimental results of Bragg. 



distributions for this airfoil at AOA values near 13.6°, as seen in figures 30-32, docs not 
indicate a collapse of the pressure peak near the leading edge nor is a trailing edge separation 
indicated. It is suspected that some error in the evaluation of the data has resulted in this low 
value. As a consequence of this experimental behavior, the ARC2D calculations were 
performed in the steady-state mode for AOA values up to 13°. 

If, however, the value of C T were desired, a time-accurate calculation at larger 

AOA would be required. This was attempted using both models and resulted in the upper part 
of the curves shown in figure 28. As seen, the Baldwin-Lomax model produces a 
of 18° AOA while the MML model yields a value of 14*. It would be interesting to determine 
which prediction produces a more accurate value. In any case, this is a first attempt 

at prediction of iced airfoil performance prior to experiment. 

The drag values shown in figure 29 indicate good agreement between calculation and 
experiment. The two turbulence models agree well except at large AOA values where the 
MML model appears to give higher values. This is due to earlier separation predicted by the 
MML model. The Cj^ values at low AOA’s are slightly underprcdicted by both calculations. 
Poor resolution of the near-wall behavior is most likely at fault. Use of a smaller grid spacing 
near the wall is suggested. It is anticipated that alteration of the C 2 value in equation (4.13) 
of the MML model may also affect the velocity gradient since this value alters the rate at 
which the mixing length approaches its limiting value. Future studies of these effects are 
suggested in order to verify this speculation. 

The type of stall that occurs for this airfoil and ice shape is indicated by examination 
of the Cp distribution and stream function contours for results past max - The pressure 
distribution shows no drop in the pressure peak near the leading edge, as seen in figure 33. 
There is however a significant bulging of the Cp values near the trailing edge. Thus, a trailing 
edge separation is indicated. In fact, examination of the stream function contours, figure 34, 
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Figure 30 Pressure coefficient distribution for NACA63^- 415 airfoil with R7 ice shape. 
Moo = 0.14, Re = 5.0x10 s , AOA = 12.6‘. Bragg [70] 
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Figure 33 Pressure coefficient distribution for NACA63^- 415 airfoil with R7 ice shape. 
Moo = 0.14, Re = 5.0x10 s , AOA = 16'. 
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plainly indicates the existance of a reverse flow region at the trailing edge. This behavior is 
evident with either turbulence model. However, the Baldwin-Lomax model suppresses this 
behavior until 18° AOA, while the MML model indicates separation at a lower AOA. 

The contours of eddy viscosity at these AOA’s indicate why this occurs. The 
Baldwin-Lomax model produces (fi t / ^i)max values of approximately 2000, as seen in figure 
35(a), with the largest values being centered inside the reverse flow region. The MML model 
produces ( fi t /fi)max values of approximately 1000, as seen in Figure 35(b), with the largest 
values being centerd at the separation region. The larger values of the Baldwin-Lomax model 
tend to suppress the development of these separation regions and results in delayed prediction 
of stall. The relaxation of eddy viscosity values aft of the separation point, as found with use 
of the MML model, allows for development of the reverse flow region. This in turn, results in 
prediction of the physically correct lower stall angle. 

NACA63^-415, G3 ice shape 

This shape was also tested in the NASAIRT wind tunnel with nominal Reynolds and 
Mach numbers of 5xl0 6 and 0.14 respectively. The data was taken for AOA values up to 
11.6’. The experimental drag values are only available for AOA values up to 7.6 , due to 
limitations in the test procedure. Comparisons between data and computations for and 
Cp are shown in figures 36 and 37 respectively with computed results for both the Baldwin- 
Lomax and MML models included. 

The lift coefficents do not agree well with the experiment for this shape. Both 
models predict lower values than experiment for AOA values from 6° to 12°. The Baldwin- 
Lomax model starts to yield values of C L higher than the experimental value at 12° AOA. 
The computed results indicate separation at this AOA. As in the case of the R7 ice shape, the 
two turbulence models yield significantly different flowfields. The Baldwin-Lomax model 
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produces a much smaller separation region than the MML model. This is again due to the 
differences in magnitude and distribution of eddy viscosity values as calculated by the alternate 
models. The pi % values determine the size of the boundary layer, which in turn influences the 
pressure distribution on the airfoil and hence the force coefficient values for these flow 
conditions. Additionally, these dissimilar turbulent viscosities produce different values of Cp 
which can result in further differences of the force coefficients. 

The computational results, using both turbulence models for AOA at Cj , are 
compared to experiment in Table 5. 



Ct 

1j max 

AOA 

Experiment 

1.2 

8.6 

Baldwin - Lomax 

1.25 

18.0 

MML Model 

0.94 

12.0 


Table 5 Calculated stall conditions using either turbulence model. 
Comparison to experimental conditions. 


As these values indicate, the stall angle is more accurately evaluated by use of the MML 
model. The value of C T is found by the Baldwin-Lomax model but this would seem to be 

Li max 

fortuitous, since it is at the wrong angle of attack. The reason for this result is somewhat 
different than for the previous case. The maximum values of /i t obtained with the MML 
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model are approximately 2500 and are distributed throughout the flowfield, as seen in figure 
38(b). The Baldwin-Lomax model, on the other hand, produces values which are close to 10 in 
some regions (i.e. near the leading and trailing edges) and close to 10,000 in other regions (i.e. 
at mid-chord locations), as seen in figure 38(a). This large variation in fi t values corresponds 
to the motion of the separation bubble over the surface. The small values are due to selection 
of Fmax at a very small y value, while the large values are due to selection of F maT at a [ 

very large y value. The MML model, on the other hand, tends to correctly concentrate // t 
values in regions of high vorticity. This results in a more realistic distribution of /i f 
throughout the flowfield. A similar fi t distribution is seen in the results of Majumdar and 
Rodi [71] for circular cylinders. f 

Interestingly, calculations using either turbulence model indicate the same mode of | 

stall behavior. This is seen in figure 39(a) which shows the stream function contours for 8° j 

AOA, that is, just prior to stall. The figure shows two separation regions on the airfoil at both ■ 

the leading and trailing edge. This corresponds to the Cp distribution obtained by Bragg, as 
shown in figure 40. As the AOA increases, the size of both these regions increases. Finally, at 
10°-12° AOA, the two recirculation regions join and stall occurs, as shown in figures 39 (b) and : 

(c). These results are from the use of the MML model. The Baldwin-Lomax model produces 
similar results, except at higher AOA values. | 

The C p distribution obtained using the MML model at 8° AOA is also shown in 
figure 40. As seen, the pressure distributions indicate that the calculated leading edge 

i 

separation bubble is smaller than the measured bubble. This in turn alters the pressure 
distribution downstream of the leading edge. Evidently this leads to a larger calculated 
separation region at the trailing edge than actually occurs. Also, the influence of the trailing 
edge separation appears to alter the lower surface pressure distribution in that region. All of 
these differences lead to a much lower value of lift than measurements indicate. 
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Figure 38 Eddy viscosity contours for the NACA63^- 415 airfoil with G3 ice shape, 
AOA = 10°, (a) Baldwin-Lomax model, (b) MML model. 
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Figure 39 Stream function contours for the N AC A 63^- 415 airfoil with G3 ice shape. 

(a) 8’ AOA, (b) 10* AOA, (c) 12' AOA 
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The inaccuracy of the Cp values near the leading edge may be produced by a lack of 
grid resolution in the separation bubble region. The complex nature of the flowfield in this 
region requires a significant number of points in both the streamwise and transverse directions. 
The measured pressure values imply a local edge velocity somewhat lower than for an attached 
flow case. The calculated values do not reflect this behavior and indicate that the local edge 
velocities are close to those for attached flow. Yet, the stream function contours in figure 39(a) 
plainly indicate that a separation region has been calculated. The measured results also imply 
that the velocity gradient is close to zero in this region. The calculated results indicate a large 
value for the velocity gradient and thus a strong reverse flow region is suspected. Comparisons 
of velocities in this region are not available for this airfoil, but they were measured for the next 

case and will be discussed further in the next section. 

Despite the differences shown in the lift values, the drag results show a remarkable 
degree of agreement. At the largest AOA’s, the Baldwin-Lomax model produces drag values 
higher than the measured values. This is due to the large values of ft, calculated by this 
model. The pressure drag is not a major factor for the Baldwin-Lomax calculations since the 
flow has not separated at the AOA values shown in the figure. The lower values of the MML 
model are a result of lower p t values and seem to reproduce the measured values very well. 
The inaccuracies in modeling the separation bubble may be the reason for disagreement 
between the MML model and experiment at these lower AOA values. At higher AOA values, 
the flow separates and the C D values are dominated by pressure forces. Experimental drag 
values were not obtained at these higher angle of attack conditions and thus no comparisons 
could be made. 

The results of this section indicate that the difficulties in modeling separated flow 
behavior may lie in the grid resolution in those regions where separation occurs, such as the 
leading edge horn, and in the turbulence modeling. More detailed information on the velocities 
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in these regions could help in diagnosing the problems in the calculations. A more detailed 
dataset was obtained in the test by Bragg [66] on a NACA 0012 airfoil. Results of the ARC2D 
calculations will be compared to Bragg’s data in the next section. 

NACA0012, G1 ice shape 

This artificial ice shape was tested by Bragg [66] at the Ohio State subsonic wind 
tunnel. His tests were conducted at a nominal Reynolds number of 1.4xl0 6 and a nominal 
Mach number of 0.12. The results of this experimental program included lift and drag data as 
well as pressure distributions and velocity profiles. The lift values were obtained by 

integration of the pressure data over the surface of the airfoil. The drag data was obtained 

from a total pressure survey made in the airfoil wake. Comparisons between the experimental 

C^and Cp values and the computed results are shown in Figs. 41 and 42, Both the Baldwin- 

Lomax and MML models were used. 

The lift results indicate the differences between the two turbulence models. The 
Baldwin- Lomax model predicts continued increase of the lift past the experimental C T 

L mar 

value. At 8° AOA, the MML model predicts unsteady vortex shedding while the Baldwin- 
Lomax model predicts a large attached recirculation bubble resulting in the higher value. 
At 10° AOA, both models predict unsteady behavior but with different and C ^ values. 
This is due to the large value of the pressure peak associated with the recirculation region 
being shed from the upper surface. The pressure peak predicted by the MML model is smaller 
and when integrated with respect to time produces more reasonable lift and drag values for the 
ice shape/airfoil combination. 

The size of the recirculation region can be altered by changing the values of C) and 
C 2 in the MML model, as shown in the study of the clean NACA0012 airfoil. This indicates 
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Figure 41 Lift coefficient vs. AOA for NACA0012 airfoil with G 1 ice shape. 
Comparison of ARC2D results using both turbulence models to 
experimental results of Bragg. 
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that correct determination of the vortex size and shedding frequency may lead to a better 
evaluation of the performance characteristics of the iced airfoil near stall conditions. 

In an effort to evaluate this effect, the C, values of 1000 and 2000 were used for 
several AOA’s near stall. The results for 8” AOA are shown in figures -13-10, which indicate 
pressure coefficients and lift histories for the two C, values. The lower C, value results in a 
larger vortex, as seen in a comparison of figures 43 and 45. The shedding frequencies are 
indicated in figures 44 and 46. These results show that the C, value of 1000 yields a time 
averaged lift coefficient of 0.46 and a value of 2000 yields a time averaged lift coefficient of 
0.47. The experimental value, taken from Bragg [66] for this AOA, was found to be 0.54. 
Thus, this term does not seem to influence the magnitude of the lift at this angle of attack. 

However, examination of the lift histories reveals some differences due to C, values. 
These plots indicate a periodic behavior of the lift. Contour plots of the stream function, 
figure 47, show the development of a large recirculation region originating at the leading edge, 
moving along the upper surface of the airfoil and eventually shedding off the trailing edge of 
the airfoil. The shedding frequency is characterized by the Strouhal number and was found to 
be 0.0100 and 0.0088 for each C, value, respectively. The work of Zaman and Mckmzie [68] 
indicates that there is a low frequency oscillation in the flow over an airfoil at conditions near 
stall. They speculate that this oscillation is due to the periodic formation and breakdown of a 
large separation bubble. The frequency of the oscillation resulted in a Strouhal number of 
0.02, which is an order of magnitude lower than the normal value associated with bluff- body 
shedding, but surprisingly close to the values mentioned above. They state, that at higher 
incidence angles they are able to produce the more conventional bluff-body shedding frequency. 
Their results suggest another area of potential computational investigation, that is, the higher 
frequency shedding at larger AOA’s. 
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Figure 44 Lift history for NACA0012 airfoil with G1 ice shape. 

Moo = 0.12, Re = 1.4xl0 6 , AOA = 8”. C, = 1000 
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Figure 46 Lift history for NACA0012 airfoil with G1 ice shape. 

Moo = 0.12, Re = 1.4xl0 6 , AOA = 8*. C x = 2000 


108 




Figure 47 Stream function contours for the NACA0012 airfoil with G1 ice shape 
AOA = 8". 
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The computational results for the iced airfoil correspond to the description given by 
Zaman and McKinzie. Thus, the Strouhal number may be used to adjust the and C 2 
values for a given airfoil geometry. The criterion is that at AOA values near stall the Strouhal 
number should be approximately 0.02 and at higher values of AOA the Strouhal number 
should be approximately 0.2. The point at which to switch from requiring a value of 0.02 to a 
value of 0.2 is not readily apparent. This is an area requiring further investigation both from 
the computational and experimental viewpoints. 

The behavior at higher AOA values was examined in order to obtain some insight 
into this problem. At 10° AOA, the Strouhal numbers resulting from use of Cj values of 1000 
and 2000 are 0.016 and 0.01T respectively. These appear to be in the range associated with the 
low frequency shedding phenomena. It is not clear, from these results, which C, value is 
preferred. The time-averaged lift values for these two cases are 0.576 and 0.460 respectively. 
The lift measured by Bragg was 0.479 at 9” AOA with the slope of the lift versus alpha curve 
being negative. Thus it would seem that for this case the Cj value of 2000 is preferred. 

These results indicate that the vortex shedding is essentially an inviscid phenomena 
which is modified by the presence of viscosity in the flowfield. Referring back to the equation 
for vorticity generation (i.e. Eq. 4.8), these results suggest that the mechanism for vorticity 
generation must be the severe pressure gradient near the horn. The resulting i circulating 
region which develops enhances the turbulence levels, which in turn leads to increased 
dissipation of the energy in that region. If the increase in /i ( is too large then the development 
of the recirculation region is impaired. This leads to the alteration in Strouhal number 
indicated above. Similarly, if the viscous dissipation is too low, then the recirculation region 
grows too rapidly and the resulting force coefficients are larger than expected, as indicated in 
the 10° AOA case. The questions then are: how well does the MML model predict stall 
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behavior and what values of Cj and C 2 are appropriate? A study of these effects is being 
planned and will serve as the basis for continuing work on this problem [72]. 

Local Results - Pre Stall Conditions 

Pressure coefficient distributions for 0°, 2°, and 4° AOA are shown in Figures 48-50. 
Comparisons between computed and experimental results indicate substantial agreement, with 
the exception of the region near the glaze ice horn. The computed results have a large pressure 
spike in this region which is not found in the experimental results. The pressure taps in the 
experiment are spaced at every one percent chord starting at the tip of the horn. Thus, it is 
unlikely that the spike is there and is not being measured. Examination of the computational 
grid indicates that the tip of the horn is represented by three grid points in a triangular 
arrangement. The fluid is thus forced to turn a corner which has an angle greater than 180°, 
which is not the case for the actual ice shape. This suggests that further grid refinement may 
be required to eliminate this pressure spike. An alternate grid code is presently being 
evaluated for this purpose, but was unavailable for use in this investigation. 

The experimental results show a flat pressure profile corresponding to the 
recirculation region. This means that the flowfield is adjusting to the concave region aft of the 
horn by filling it with low-velocity recirculating fluid. The shear layer which lies on top of this 
stagnant region is thus flowing over a virtual surface which is more aerodynamic than the 
actual iced airfoil geometry. The computer code captures this recirculation region but does not 
adjust the pressure field appropriately. 

The code does allow variation of the pressure in the direction normal to the surface, 
however, it does not seem to capture this feature of the flow. This is indicated in figure 51, 
which shows the static pressure contours in the recirculation region. The low pressure levels 
found near the horn do not extend into the recirculation region as the experimental results 
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AOA = 2*. (a - upper surface, b - lower surface) 
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Figure 51 


Static pressure contours in the recirculation region behind the horn 

for 4° AOA. 
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would suggest. The experimental results do imply that the pressure gradient in the ^-direction 
may be as important as the ^-direction gradient. Hence the viscous terms in the ^-direction, 
especially the v f terms, may also be of some importance. These terms are neglected in the 
thin-layer form of the Navier-Stokes equations. The ARC2D code has the capability to include 
the explicit portion of these terms. When this option is employed, there is no significant 
difference to the solution. This suggests again that the grid resolution in this region is 
insufficient. It is expected that an alternate grid will resolve this problem along with the 
pressure spike at the tip of the horn. 

The velocity profiles also suggest a need for greater grid resolution in the recirculation 
region. The velocity profiles obtained using ARC2D are compared to experimental values in 
figures 52-54. These figures show the velocity profiles within the recirculation region for AOA 
values of 0”, 2°, and 4° respectively. Notice the differences in height of the zero-velocity point 
and in location of reattachment. The height of the reverse flow region, defined as the distance 
from the surface to the zero-velocity line, can be as large as 2-3 percent chord. The grid 
resolution in this region is not as fine as it is near the surface. This is true for resolution in 
both the r) and £ directions. The results presented in figures 52-54 were obtained using the 
Baldwin-Lomax model. The turbulence model selection does not seem to play as important a 
role as either the transition location or the grid resolution. Further investigation of the 
turbulence model role in the development of the attached bubble is required. 

The free shear layer, extending from the separation point to the reattachment point, 
should have the same grid resolution as the boundary layer. The fact that it docs not means 
that certain aspects of the flow physics are not being modeled correctly. This is also reflected 
in figures 52-54, by examining the velocity gradients near the upper edge of the reverse flow 
regions. The code tends to underpredict the velocity gradient along with the distance of the 
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Figure 52 Velocity profiles in the recirculation region aft of the horn 
ARC2D, o Experiment, AOA = 0* 
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Figure 53 Velocity profiles in the recirculation region aft of the horn 1 

ARC2D, o Experiment, AOA = 2* * 
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free shear layer from the surface. This inability to capture the shear layer behavior is 
consistent with the overprediction of the pressure spike mentioned previously. 

Correct representation of the flow in these recirculation regions also requires 
specification of the transition region. Presently, the ARC2D code has a very rudimentary 
single-location transition specification. As described by Mehta [73], ’the process of transition 
from laminar to turbulent flow is not a well-defined problem because of the sensitivity to 
poorly defined initiating disturbances’. He goes on to state that the main requirements for 
numerical solution to the compressible Navier-Stokes equations in regard to the transition 
process are ’that (1) the discretization errors do not contaminate physical phenomena such as 
instabilities, that is, the finest scales in the transition process are adequately resolved in space 
and time; and (2) the introduction of artificial boundaries owing to the limited size of the 
computational domain does not interfere with the physical upstream influence, the ellipticity of 
the Navier-Stokes equations.’ 

From these comments, it is apparent that appropriate representation of transition is 
necessary for correct modeling of the recirculation region aft of the glaze horns. Since the code 
presently does not have a sophisticated transition model, the effect of transition specification 
was examined by simply altering the location and noting the changes in the velocity profiles. 
This was done for the 0° AOA case by moving the transition location from the tip of the horn, 
as was the case in figures 52-54, to a point approximately in the center of the separation 
bubble. The results are shown in figure 55. 

Moving the transition location further aft has two effects on the separation bubble. 
The modeling of the shear layer velocity gradient is improved and the magnitudes of the 
reverse flow velocities are overpredicted. The improvement to the shear layer velocity gradient 
can be attributed to a more realistic representation of the dissipation in that region. The use 
of the downstream transition location more accurately represents the free shear behavior, as 
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Figure 55 Velocity profiles in the recirculation region aft of the horn 
ARC2D, o Experiment, AOA = 0" 
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described by Mehta. The overprediction of the reverse flow velocities is due to the approach 
used by ARC2D in defining laminar and turbulent regions. 

The code first calculates turbulent viscosities for all constant-^ lines from the leading 
edge to the trailing edge. The transition points on the upper and lower surfaces are then 
identified and all turbulent viscosity values from the leading edge to these transition points arc 
reset to zero. If the transition location is in the center of a separation bubble then turbulent 
fluid re-laminarizes as it flow toward the forward portion of the bubble in the reverse flow 
region. This decrease in dissipation leads to an inappropriate enhancement of the velocities in 
that region. 

The behavior described by Simpson et al. [63] for a two-dimensional turbulent 
separated flow seems more realistic. He indicates that the mean backflow in the detachment 
region of the separation bubble is a result of incursions of turbulent fluid from the overlying 
shear layer. The use of a laminar flow region in this part of the bubble is inappropriate. Also, 
the use of a turbulence model based on the mean flow velocity profiles seems precluded by this 
description. Thus, at this time no readily available method will adequately describe the 
detachment region. The lack of such a model requires the selection of some transition location 
which produces the most acceptable results from the standpoint of performance evaluation. It 
seems that a transition location halfway through the separation region produces a better shear 
layer evaluation and reattachment point prediction. Thus, selection of this point for transition 
location is recommended at this time. Future work could be directed at improvement of the 
model with respect to the type of intermittant behavior described by Simpson. 

Local Results - Post Stall Behavior 

The ice accretion geometry and angle of attack play a critical role in altering the 
maximum lift and inducing the onset of stall. The results of flow visualization studies 
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indicate, for the NACA0012 airfoil and G1 ice shape, that the separation bubble near the horns 
grows with increasing incidence until at approximately 6* — T AOA the entire upper surface 
exhibits alternating forward and reverse flow. Results from the computations also indicate a 
significant change at 6° — 7° AOA. As mentioned earlier, the MML model predicts unsteady 
flow at 6° AOA while the Baldwin-Lomax model predicts the onset of unsteady flow at 10° 
AOA. Therefore, the MML model will be used exclusively for evaluation of the post stall 
behavior described in this section. 

At lower AOA’s, the separation bubble size is a function of AOA but it exhibits an 
essentially steady size and shape. At 6° AOA, the bubble grows to encompass the entire 
surface and continues growing in strength until it is swept from the surface by the main flow. 
This flow is periodic in nature with new bubbles being generated at the leading edge to replace 
those being shed at the trailing edge. This shedding process is present for a large range of 
AOAs above stall. This sequence is captured in figures 56-74, for a 10° AOA case. The lift 
history for this process is shown in figure 75. 

These figures show the development of the bubble at various stages of growth and 
subsequent shedding over several cycles. The figures shown are the stream function and equi- 
vorticity contours at selected time points in the computation. .Initially, the separated flow 
region encompasses the entire upper surface of the airfoil, as seen in figure 56. The equi- 
vorticity contours, on the other hand appear more concentrated near the horn. This is the 
point in the process just after one shedding cycle and just prior to the next cycle. At this 
point, the lift of the airfoil is just past its maximum. The zero stream function contour has 
just separated from the surface and the lift is starting to decrease. A small region of counter- 
clockwise (i.e. positive) rotating fluid is seen at the trailing edge. The high lift value is due to 
the large amount of circulation within the larger, negative vortex which is still present on the 
upper surface. 
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Figure 57 Vortex development on the NACA0012 airfoil with G1 ice shape; AOA=8”, 
t=t 0 +300At; (a) stream function contours, (b) equi-vorticity contours 
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Figure 58 Vortex development on the NACA0012 airfoil with G1 ice shape, AOA 8 , 
t=t 0 + 600 At ; (a) stream function contours, (b) cqui-vorticity contours 
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(b) 


Figure 59 Vortex development on the NACA0012 airfoil with Gl ice shape; AOA=8°, 
t=t 0 + 900 At ; (a) stream function contours, (b) equi-vorticity contours 
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Figure 60 Vortex development on the NACA0012 airfoil with G1 ice shape; AOA 8 , 
t=t o +1200At; (a)stream function contours, (b)cqui-vorticity contours 
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Figure 61 Vortex development on the NACA0012 airfoil with G1 ice shape; AOA=8\ 
t=t 0 + 1500At ; (a) stream function contours, (b) equi-vorticity contours 
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Figure 63 Vortex development on the NACA0012 airfoil with G1 ice shape; AOA=8°, 
t=t o +1700At; (a) stream function contours, (b) equi-vorticity contours 
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(b) 


Figure 64 Vortex development on the NACA0012 airfoil with G1 ice shape; AOA-8”, 
t=t o +1800At; (a) stream function contours, (b) equi-vorticity contours 
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Figure 65 Vortex development on the NACA0012 airfoil with G1 ice shape; A0A=8 o , 
t=t 0 + 1900At; (a) stream function contours, (b) equi-vorticity contours 
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Figure 66 Vortex development on the NACA0012 airfoil with G1 ice shape; AOA=8°, | 

t=t o -f2000At; (a) stream function contours, (b) equi-vorticity contours ; 
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Figure 67 Vortex development on the NACA0012 airfoil with G1 ice shape; AOA=8\ 
t— t o -|-2100At; (a) stream function contours, (b) equi-vorticity contours 
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Figure 68 Vortex development on the NACA0012 airfoil with G1 icc shape; A0A=8 , 
t=t o +2200At; (a) stream function contours, (b) equi-vorticity contours 
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Figure 69 Vortex development on the NACA0012 airfoil with G1 ice shape; AOA=8°, 
t=t 0 -f 2300At; (a) stream function contours, (b) equi-vorticity contours 
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Figure 70 Vortex development on the NACA0012 airfoil with G1 ice shape; AOA=8\ 
t=t 0 -f 2400 At; (a) stream function contours, (b) cqui-vorticity contours 
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Figure 71 Vortex development on the NACA0012 airfoil with G1 ice shape; AOA=8°, 
t=t o +2500At; (a) stream function contours, (b) equi-vorticity contours 
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(b) 


Figure 72 Vortex development and shedding on the NACA0012 airfoil with G1 ice shape; 

t=t 0 + 2600 At ; (a) stream function contours, (b) equi-vorticity contours 


140 
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Figure 73 Vortex development and shedding on the NACA0012 airfoil with G1 ice shape; 

t=t 0 + 2700At ; (a) stream function contours, (b) equi-vorticity contours 
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Figure 74 Vortex development on the NACA0012 airfoil with G1 ice shape; AOA=8”, 
t=t o +2800At; (a) stream function contours, (b) cqui-vorticity contours 
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The next figure in the sequence has two significant features. The leading edge region 
exhibits a pinching off of the separation bubble and a gradual movement of the vortex center 
downstream and away from the surface. This coincides with a distinct shedding of vorticity 
from the surface toward the free stream. Also, the small bubble near the trailing edge has 
merged with the larger bubble and the resulting bubble has reattached to the surface. This 
second phenomenon is a residual adjustment to the previously shed vortex. 

The bifurcation of the leading edge separation region can be attributed to the 
convection of vorticity away from the horn, where it is initially developed. As the process 
starts, a vortex is present just past the ice shape region. This vortex moves away, as time 
progresses, due to entrainment by the shear layer. It is replaced by vorticity which is 
constantly being generated at the horn itself. This last point is indicated by the constant 
vorticity level found at the horn, even as vorticity is being converted away. This is seen in 
figures 76-85, which show equi-vorticity contours near the horn. Thus, as one separation 
bubble moves along the airfoil surface another is being created at the horn to lake its place. 
The timing of these two processes determines the amplitude and frequency of the lift and drag 
oscillations. 

The reattachment of the bubble at the trailing edge keeps the separation bubble on 
the airfoil. This prevents the lift loss from being more severe. The vortex associated with this 
bubble can be clearly seen in figure 57. The drop in lift of the airfoil can be associated with 
the movement of this vortex away from the surface and into the wake, shown in figures 57-59. 
These figures show the growth of the negative and positive separation bubbles at the leading 
and trailing edges, respectively. These bubbles grow in strength and force the larger bubble to 
first separate from the surface (fig. 58) and then to flow into the wake (fig. 59). This vortex 
sustains lower pressures on the upper surface and when it leaves the surface these pressures 
increase. This collapse of the pressure peak produces the drop in lift. 
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Figure 76 Equi-vorticity contours near the leading edge of the NACA0012 airfoil 
with G1 ice shape; AOA=8°, t=t 0 + 1500 At 


145 




Figure 78 Equi-vorticity contours near the leading edge of the NACA0012 airfoil 
with G1 ice shape; A0A=8°, t=t 0 4- 1700At 
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Figure 79 Equi-vorticity contours near the leading edge of the NACA0012 airfoil 
with G1 ice shape; AOA— 8°, t=t 0 -f 1800At 
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Figure 80 Equi-vorticity contours near the leading edge of the NACA0012 airfoil 
with G1 ice shape; AOA=8°, t~t 0 + 1900 At 
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Figure 81 Equi-vorticity contours near the leading edge of the NACA0012 airfoil 
with G1 ice shape; AOA=8\ t~t 0 + 2000 At 
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Figure 82 Equi-vorticity contours near the leading edge of the NACA0012 airfoil 
with G1 ice shape; AOA=8°, t=t 0 + 2100 At 
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: Figure 83 Equi-vorticity contours near the leading edge of the NACA0012 airfoil 

i with G1 ice shape; A0A=8 , t=t 0 + 2200At 

i t 
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Figure 84 Equi-vorticity contours near the leading edge of the NACA0012 airfoil 
with G1 ice shape; A0A=8°, t=t 0 + 2300At 
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Figure 85 Equi-vorticity contours near the leading edge of the NACA0012 airfoil 
with G1 ice shape; A0A=8°, t=t 0 + 2400 At 
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There are two other interesting features of this shedding process. The vorticity on 
the lower surface remains undisturbed by the activity on the upper surface. The vorticity 
levels here are established by the airfoil/ice shape geometry and the angle of attack. 
Additionally, the leading edge stagnation point remains at the same location throughout the 
process. The stagnation region appears to be well isolated from the unsteady flow by the horn 
of the ice shape. This last point has significant consequences for the modeling of the ice 
accretion process itself. It appears that extreme accuracy in modeling of the flowfield 
downstream of the ice shape is not necessary for correct representation of the incoming flow. 
Some accounting for separation may be necessary, however, for establishing the correct limits 
of water droplet impingement. Certainly, further examination of the flow near the separation 
region should be undertaken in order to confirm this assertion. 

Figure 60 indicates that the vorticity generated at the horn has again grown to 
encompass the entire upper surface of the airfoil. This occurs as a result of the movement of 
the previous bubble off the airfoil surface. The flowfield pattern is now similar to that of 
figure 56, however, there are significant differences which result in a slightly altered shedding 
sequence for this bubble. In this case, the bulge in the separation bubble is centered more 
toward the trailing edge. This would suggest a shorter shedding period, as the bubble is 
already further along the surface. Indeed, figures 61-65 show that half the bubble has moved 
into the wake before the positive bubble develops at the trailing edge. This produces a 
somewhat different shedding behavior and consequent changes in lift and drag. 

During this shedding period, the large negative bubble is not shed completely but is 
pinched off by the positive bubble. Figures 65-67 show how the negative bubble is essentially 
cut in two with some fluid moving into the wake and some being forced back onto the airfoil 
surface. As a result, the circulation on the surface is restored as indicated by the reattachment 
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of the bubble seen in figure 68. This causes the decrease in lift to stop at a Cj value of 
approximately 0.3 as compared to the minimum value of -0.1 from the previous shedding 
cycle. 

Examination of figures 76-83 indicates that during this time a new negative vortex 
has developed at the horn. This vortex grows in size and strength until, during the interval 
from t = t 0 ~f 1800 At until t = t o +2200At, the vortex separates from the horn and moves 

F 

further back along the airfoil surface. This event corresponds to the detachment of the 
separation bubble from the horn. The vorticity level at the horn itself remains essentially 

f 

constant, apparently independent of the shedding process occurring downstream. The vorticity \ 

level at the horn appears to be established by the horn geometry, the angle of attack, and the I 

freestream conditions. Consequently, modeling the flow in this region is extremely important 

for establishing the correct size, strength, and timing of the vortex shedding process. How 

much detail of the horn geometry is required, what grid size and spacing is appropriate, and 

what type of flow modeling is appropriate are all pertinent questions when trying to capture : 

the details of this process. 

Figures 68-74 show the movement of the separation bubble along the airfoil surface 

and its subsequent shedding into the wake. This case is more like the first shedding event, in 

that the separation bubble is almost completely off the surface by the time the positive bubble 

develops. It is interesting to note the changes occuring to the separation bubble, as it moves 

along the airfoil surface. More fluid appears to enter the bubble as indicated by the larger 

number of stream function contours. At the same time, the vorticity levels do not seem to 

change considerably, thus indicating no significant change in the amount of circulation | 

| 

associated with this bubble. The lift increase must, therefore, be due to the vorticity being 
generated at the horn. When the bubble finally bursts (fig. 72), the drop in circulation is 

■ 

much greater than the amount being generated at the horn and hence the lift decreases rapidly. | 
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The positive vorticity which develops at the trailing edge is seen to remain on the 
airfoil upper surface (i.e. figure 74). This further reduces the lift such that, at t o -}-2800At, it is 
negative. This region of positive vorticity eventually merges with the leading edge vortex, as 
was the case in figures 59 and 60. This merging of the two regions and the continuing creation 
of negative vorticity at the leading edge leads to an eventual restoration of the lift, until the 
entire cycle repeats. 

The turbulence model can have a significant impact on the development of the 
separation bubbles associated with this shedding process. As seen in the previous chapter, the 
Baldwin-Lomax model has some difficulties with separated flow due to the specification of the 
F max parameter. The behavior of the MML model during the vortex shedding is seen in 
figures 86-91. The regions of highest turbulent dissipation correspond to locations of 
separation and reattachment. This insures that velocity levels in these regions remain low, 
consistent with expectations. These regions move along the surface, following the motion of 
the separation bubbles. As the bubble moves along the surface, the turbulence level decreases 
reflecting the lower mean flow velocities within these regions. 

The other interesting feature of the turbulence model behavior during this sequence is 
the development of the leading edge region. The eddy viscosity levels just aft of the 
upper-surface horn remain consistently high throughout the shedding cycle. Thus, as vorticity 
is constantly generated at the horn the turbulence level remains correspondingly high. These 
/x, levels can in turn influence the growth rate of the vorticity and some type of feedback 
mechanism is thus established. Accurate prediction of turbulence levels near the horn is 
therefore essential to the correct prediction of the vortex shedding time scales and hence of the 
integrated force coefficient values for this post-stall behavior. 

Finally, examination of the lower-surface indicates steady flow behavior with 
associated constant eddy viscosity values. The vorticity generation rate at this horn must be 
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0 



(b) 


Figure 86 Turbulence model behavior during vortex shedding process; AOA=8", 

t=t o +2400At; (a) stream function contours, (b) eddy-viscosity contours 


158 



(b) 


Figure 87 Turbulence model behavior during vortex shedding process; AO A— 8°, 

t=t o T2500At; (a) stream function contours, (b) eddy- viscosity contours 
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(b) 


Figure 88 Turbulence model behavior during vortex shedding process; AO A =8", 

t=t o +2600At; (a) stream function contours, (b) eddy- viscosity contours 
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(b) 


Figure 90 Turbulence model behavior during vortex shedding process; AOA=8", 

t=t o -f2800At; (a) stream function contours, (b) eddy- viscosity contours 
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(b) 


Figure 91 Turbulence model behavior during vortex shedding process; AOA=:8°, 

t=t 0 T 2900At; (a) stream function contours, (b) eddy-viscosity contours 
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lower than for the upper-surface and can be easily convected along the airfoil. This lower 
generation rate precludes the development of large regions of vorticity and hence the lower- 
surface does not experience the unsteady behavior seen on the upper-surface. The turbulence 
levels on this surface are substantially lower than those of the upper-surface. Tins is consistent 
with the near laminar behavior on the pressure surface of a clean airfoil. The favorable 
pressure gradient on the lower-surface does not promote the development of a large amount of 
vorticity away from the surface. This lower vorticity level in turn produces a lower eddy 
viscosity level. In fact, the turbulence on this surface appears to develop at the lower- surface 
horn and is then distributed along the airfoil. 

As stated previously, the difference in time scales of vorticity generation and 
convection lead to alteration in frequency and amplitude of the force coefficient fluctuation. 
Since the eddy viscosity level can affect both of these time scales, the turbulence model can 
influence the temporal development of the lift, drag, and pitching moment. Correct 
representation of the turbulence level is therefore essential to the modeling of the post-stall 
behavior of an airfoil, either clean or with leading edge ice. While the MML model may not 
provide completely accurate predictions of turbulence level for a stalled airfoil, it can provide 
information on the appropriate turbulence levels and distributions for modeling of such 
behavior. 

Further investigations will be required to determine the sensitivity of t lie MML model 
to different geometries and grid resolutions. Airfoil geometries can result in several types of 
stall behavior; leading-edge stall, trailing-edge stall, and combined leading and trailing edge 
stall. The ability of codes such as ARC2D to calculate stall behavior may also depend on 
which stall mechanism occurs. For the iced airfoil, the mechanism is more apparent than in 
most clean airfoil stalls. The iced airfoil case may result in a greater degree of two 
dimensionality due to the highly structured behavior near the horn. If this is the case, the iced 
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airfoil computations may be a truer representation of the actual physics. This is discussed, for 
the clean airfoil, by Zaman and McKinzie [68]. 

Grid sensitivity studies must be performed in order to determine requirements for 
capturing of the vortex shedding. Preliminary results have indicated a dependence of 
integrated force coefficients on the spacing near the surface and around the horn. It is 
necessary to determine the magnitude of this dependence and how fine a grid is required. 
Methods are presently being developed for adapting grids to the characteristics of the particular 
flowfield being studied. These may be required for unsteady flows such as those being studied 
here. Also, unstructured grids have the potential to provide a better representation of irregular 
geometries, such as the iced airfoil. These and other developments show great promise in 
providing better methods for iced airfoil analysis in the future. 
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CHAPTER 6 


SUMMARY 


Evaluation of the aerodynamic performance changes to an airfoil due to leading edge 
ice accretions has been performed using a 2D unsteady Navier-Stokcs computer code. 1 his 
code, along with appropriate grid generator and turbulence model, has provided considerable 
information on the aerodynamics of these highly complex geometries. Comparisons to 
experimental data have provided insight into the capabilities and limitations of this modeling 
method. Calculations indicate a high degree of confidence in the results for integrated force 
coefficients at pre-stall conditions. Discrepencies were noted in the separation bubble at pre- 
stall conditions and in the unsteady flowfield at post-stall conditions. A new turbulence model 
formulation was employed in order to address these difficulties and met with some degtee of 
success. Further use of this turbulence model along with finer grids may help in evaluation <>r 
these conditions. 

The MML turbulence model employed for this study is a zero-equation, eddy 
viscosity model with a Modified Mixing Length formulation which does not incorporate the 1 - 
function of the Baldwin-Lomax model and does not require the calculation of the displacement 
thickness. In this model, the mixing length is based on the local value of y* and the velocity 
scale is based on the mixing length and the local value of the vorticity. In this way, the 
feedback existing between the wall shear and the turbulence level is not dependent on I In- 
distribution of the vorticity in the far field but rather on its value near the wall. I his pie\ents 
large scale structures, which have essentially inviscid behavior, from unduly influencing the 
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turbulence level throughout the flowfield. In this model, the level of the turbulence is 
established by the fluid behavior near the wall, where vorticity generation takes place, and its 
distribution throughout the fluid is determined by the vorticity in the far field. This allows 
the development of the large vortices seen in the post-stall conditions examined and yet docs 
not ovcr-predict their size due to unrealistically low turbulent viscosity levels. 

The use of the MML model was examined by evaluation of flow behavior on a clean 
airfoil. Calculations were performed for a clean NACA 0012 airfoil under attached and 
separated flow conditions. The attached flow results indicated that the Baldwin-Lomax model 
resulted in a more accurate representation of the drag. This was the result of better 
representation of near-wall behavior. The lift values were quite accurate for both models and 
did not have any impact on evaluation of their relative performance. The higher angle of 
attack calculations indicated qualitative differences between the two models. The MML model 
calculation indicated an unsteady flowfield with periodic vortex shedding. The Baldwin-Lomax 
model, on the other hand, indicated a steady, attached flow with a value much higher than 
experiment. Additionally, a laminar flow calculation was performed and the unsteady flow 
behavior was again observed, albeit with a different shedding frequency. The MML model and 
laminar flow calculations had shedding frequencies both of which have been observed in 
experiment. Which one is appropriate, for this airfoil and under these flow conditions, will 
require further investigation. 

The ability of the code, using either turbulence model, to predict the force coefficients 
for the iced airfoils is very good at angles of attack below stall. A rime condition and two 
glaze conditions were evaluated, with lift and drag values agreeing very well with experiment.. 
Pressure distributions also indicated good agreement with the exception of the region in the 
immediate vicinity of the separation bubble aft of the glaze ice horn. Results in the separation 
bubble indicated an overprediction of the pressure peak in the forward part of the bubble and a 
lack of the constant pressure profile through the aft section of the bubble. This inablility to 
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capture the correct, pressure profile resulted in poor representation of the velocity distribution 
within the bubble. Both turbulence models had difficulties in resolving this region and in fact 
the IBL approach of Cebeci [l] has also indicated similar difficulties. 

The role of transition location was examined and found to significantly alter the 
velocity profiles. Setting transition at the midpoint of the bubble resulted in better predict ion 
of the bubble height and reattachment location but also resulted in overprediction of the 
reverse flow velocities in the forward part of the bubble. Such limitations in the present 
implementation of the turbulence models prevented further investigation along these lines and 
further work is required. The role of grid spacing in modeling of the separation bubble is also 
being investigated. It is quite likely that the resolution in the shear layer, above the separation 
bubble, is not sufficient to capture the details of that flow, which could play a major role in 
prediction of reattaclirnent location. Also, the streamwise resolution may not be sufficient to 
resolve the details of the separation and reattachment processes. This is presently being 
addressed by implementation of an unstructured grid scheme which would be more easily 
adapted to the irregular geometry of the iced airfoil. 

At flow conditions near stall, the ability of the code to correctly predict lift and drag 
values is not as robust. There can be considerable differences in prediction of stall conditions 
and the AO A at which stall occurs. Several factors influence these calculations; grid spacing 
near the ice shape, transition location, and turbulence modeling all play a major role. The 
results indicate that the MML model appears to predict the Cj max values and stall angles 
better than the Baldwin-Lomax model. This is due to the lower /q values obtained with the 
MML model. These lower eddy viscosity values allow development of the vortex shedding 
process and yet alter the shedding frequency from the laminar-flow values. The frequencies 
obtained using the MML model correspond to the low frequency shedding observed by Znnian 
and McKinzie in experiments on another airfoil shape. Further examination of this unsteady 
behavior is presently underway in a joint computational and experimental st udy [72]. 



In summary, the approach of examining the performance degradation of an iced 
airfoil using a Navier-Stokes code has been shown to be feasible. Reasonably accurate results 
have been obtained for pre-stall conditions. Determination of the angle of attack of stall can 
be predicted using the MML turbulence model. The accurate calculation of post-stall behavior 
requires further development. The MML turbulence model was developed and tested for 
calculation of this post-stall behavior with positive results. Further refinement of this 
modeling procedure is required and recommendations for continued activities in this area arc 
outlined below. 

6.1 Future Activities 

The need for accurate evaluation of the aerodynamics associated with iced airfoils is 
twofold. A representation of the velocity field surrounding the iced airfoil is important for use 
in conjunction with ice accretion prediction codes. The velocity field can impact particle 
trajectories and heat transfer calculations, which in turn affect the prediction of ice shape 
development. The use of a code such as this also allows the prediction of performance 
degradation and the onset of premature stall. These requirements serve as justification for 
futher research into the development of the correct methods for predicting the complex 
aerodynamics associated with iced airfoils. The present study has shown areas requiring 
further investigation which shall now be summarized. 

i) A grid refinement study is required to determine the sensitivity of force 
coefficient calculations to the number and distribution of grid nodes. This 
study should examine the attached separation bubble at low angles of attack 
and the region of vortex shedding activity at angles of attack past stall. The 
use of an unstructured and adaptive grid method may be appropriate. 
Additionally, grid refinement near the ice shape itself may improve pressure 
distribution and heat transfer calculations in these regions. 
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ii) The effects of transition location and the modeling of the transition process 
itself should be examined with regard to modeling of the attached separation 
bubble. Some method of allowing the code to determine the transition location 
should be examined. The method of Michel [75], presently used in many codes, 
was developed for attached flow transition and may require modification for the 
shear layer associated with the separation bubble. 

iii) A reformulation of the turbulence model to account for the low flow conditions 
in the separation bubble and still provide for turbulent flow in the overlaying 
shear layer should be developed. This would allow for transition to tmhuleul 
flow in the shear layer without resulting in a relaminarization of the reverse 
flow region in the bubble. This should be coupled to experimental 
investigations of the turbulence characteristics within the bubble which an 
presently underway [76]. 

iv) Further examination of the unsteady behavior predicted at high angles of attack 
is required utilizing the results of the investigations just described. 
Comparisons to experiment, especially the unsteady components such as 
shedding frequencies and amplitudes, can be used to provide insight on the 
applicability of the grid, of the turbulence modeling, and of the flowfield code 
itself. Such information can be used to further refine the implementation of 
these components for the iced airfoil geometry. 

v) Implementation of a method for representing surface roughness will be required 
for investigation of its relative importance. Possibilities include an adjustment 
to the turbulence model to account for equivalent sand grain roughness or 
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implementation of a discrete roughness element model. The former requires a 
method of representing icing roughness by equivalent sand grain roughness and 
the latter requires the characterization of icing roughness heights, sizes, and 
spacing. Both approaches would require a correlation of the roughness to the 
environmental conditions prevalent during the accretion process. 

Future analysis activities will require the evaluation of icing effects on multi- 
element airfoils, swept wings, stability and control parameters, and total 
aircraft configurations. This means that an ability to predict 3D effects will be 
important. Thus, the extension of present methods to 3D codes must be 
pursued while acknowledging the considerable amount of work still required in 
less complicated 2D analysis. The use of a zero-equation turbulence model will 
be all the more important for 3D analysis due to its simplicity and low 
computational overhead. 
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Appendix 1 


Nomenclature 


A = <9E/<9Q 

A+ = van Driest damping constant = 26 

AOA = angle of attack 

a = speed of sound 


B 

C 


Cl = 
C 2 = 
C cp = 



C wk = 

D = 

D_ 

Dt 

E 

E = 

Ev 

Ev 


8F/dQ 

airfoil chord length 

MML turbulence model constant; controls mixing length saturation level 
MML turbulence model constant; controls blending region curvature 
Baldwin-Lomax turbulence model constant = 1.6 
coefficient of drag = D/(0.5/?U oo 2 s) 
friction coefficient = ?" w /(0.5/)U o c 2 ) 

Klebanoff intermittancy factor constant = 0.3 
coefficient of lift = L/(0.5/?U oo 2 s) 
pressure coefficient = (P-P oo )/(0.5pU oo 2 ) 

Baldwin-Lomax turbulence model outer region constant = 0.25 
drag force 

substantial derivative 

convection terms in Navier-Stokes equations; Cartesian coordinates (Eq. 3.14) 
convection terms in Navier-Stokes equations; curvilinear coordinates (Eq. 3.22) 
viscous terms in Navier-Stokes equations; Cartesian coordinates (Eq. 3.15) 
viscous terms in Navier-Stokes equations; curvilinear coordinates (Eq. 3.23) 
total energy 
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viscous term in Ev matrix of energy equat ion 


e 4 

F = convection terms in Navier-Stokes equations; Cartesian coordinates (Eq. 3. M) 

F = convection terms in Navier-Stokes equations; curvilinear coordinates (E<|. 3.22) 

F(y) = y |w|[l - exp(-y + /A + )] in Baldwin- Lomax turbulence model 

F = Klebanoff interim ttancy factor 

Fv = viscous terms in Navier-Stokes equations; Cartesian coordinates (Eq. 3.15) 

Fv = viscous terms in Navier-Stokes equations; curvilinear coordinates (Eq. 3.23) 

F wflfce = length x velocity scale in outer region of Baldwin-Lomax turbulence model 

f = vortex shedding frequency 

f 4 = viscous term in Fv matrix of energy equation 

h = At 

I = identity matrix 

J = Jacobian of the coordinate transformation 

K = Clauser constant = 0.0168 

L = Lift force 

X = characteristic length 

i — turbulence model length scale 

M = dS/dQ 

P = pressure 

P = constant in Poisson equation used in GRAPE code 

Pr = Prandtl number 

Q = independent variables in Navier-Stokes equations; Cartesian coordinates 

Q = independent variables in Navier-Stokes equations; curvilinear coordinates 

Q = constant in Poisson equation used in GRAPE code 

Re = Reynolds number 

S = viscous terms in thin-layer form of Navier-Stokes equations (Eq. 3.33) 
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St = 


t 


u 

u 


Vdiff = 


u 


oo 


u 


V 


V 


V 

X 


y 



a 


7 

6 


1 

e 


K 


Strouhal number 
time 

contravarient yelocity in curvilinear coordinate system 
velocity vector 

term in F wafce parameter of Baldwin-Lomax turbulence model (Eq. 3.45) 

freestream velocity 

velocity in x-direction 

wall shear velocity (Eq. 4.6) 

contravarient velocity in curvilinear coordinate system 

velocity in y-direction 

velocity scale in turbulence models 

Cartesian coordinate 

Cartesian coordinate 

boundary layer coordinate; Baldwin-Lomax - Eq. 3.40; MML - Eq. 4.5 

wall shear length scale (Eq. 4.10) 

angle of attack 

ratio of specific heats = 1.4 

boundary layer thickness 

curvilinear coordinate (nominally normal to body surface) 
constant in implicit time-differencing scheme 
von Karman constant = 0.4 


A = -§/* 

H — viscosity 

H t = turbulent viscosity 

v — kinematic viscosity; \ij p 

£ = curvilinear coordinate (nominally in streamwise direction) 
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density 

shear stress; time in curvilinear coordinate system 
constant in implicit time-differencing scheme 
vorticity vector 
vorticity 

Subscripts 

i = grid index in ^-direction 

t = differentiation with respect to t 

w = at the wall 

x = differentiation with respect to x 

y = differentiation with respect to y 

r] = differentiation with respect to rj 

£ = differentiation with respect to £ 

r = differentiation with respect to r 

oo = freestream conditions 

Superscripts 

n = iteration number 
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Appendix 2 


MML Turbulence Model Code Listing 

The subroutine listed on the following pages calculates the turbulent viscosity level 
for the flowfield obtained from the ARC2D calculation on a user-specified two-dimensional grid 
system. The subroutine uses the velocities and grid coordinates to determine length and 
velocity scales according to the equations described in Chapter 4. These length and velocity 
scales are then used to find [i v which is then used in the subsequent iteration of the velocity 
calculation procedure. The pertinant variable names are described below in order to clarify the 
relationship of the code to the actual MML model equations. 

APLUS = Van Driest damping factor; 26 

C2B = Temperature ratio for Sutherland law = 

^ oo 

C2BP = 1. + C2B 

DELTA(J) = Array of YSTAR values on the airfoil surface 
FMUtmp(K) = Viscosity value at a K location 
FI = Cj 
F2 = C 2 

GAMMA = Ratio of specific heats; j 
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J = Index in ^-direction 


JMAX = Maximum J value in grid 

JTAIL1 = First J grid line on airfoil surface 

JTAIL2 = Last J grid line on airfoil surface 

JTRANLO = J grid line location for transition on lower surface 

JTRANUP = J grid line location for transition on upper surface 

K = Index in ^-direction 

KMAX = Maximum K value in grid 

MXLNGTH = Mixing length; t 

NUMITER = Iteration number 

PRESS = Pressure; GAMI*(Q(4) - 0.5*(Q(2)**2+Q(3)**2)/Q(1)) 

Q(J,K,M) = Q; Conservative variables: M=1 =>p 

M=2 =>pu 
M=3 =>pv 
M— 4 =>e 
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RA = 1/y 


RE = Reynolds number 


SNOR = Distance normal to surface 


TAU = Vorticity at the surface; | 


TINF = 



TURMU = \i t 


VK = k 


VORT = Vorticity at a grid location 


WMU = /i w 


X = x location of grid point 


XY(J,K,M) = Metrics of transformation: M 

M = 
M 
M 


1 =>£x 

2 =>(,y 

3 =>r) x 
A 
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XYJ(J,K) = Jacobian 


YSTAR = y* 


YPLUS = y + 
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q**************************** TURBULENT VISCOSITY ************************* 

0**********************************************%****************************** 

SUBROUTINE MUTUR(JDIM,KDIM,Q, PRESS, VORT,TURMU,X,Y,XY,XYJ) 
PARAMETER(MAXJ=260) 


COMMON/BASE/ 


lJMAX, 

KMAX, 

JM, 

KM, 

JBEGIN, JEND, 

1 KBEGIN, 

KEND, 

jplus(999), 

jminus(999), jlow, jup, 

1 KLOW, 

KUP, 

JMAXOLD, 

PERIODIC, 

NP, DT, 

1 FSMACH, 

ALPHA, 

GAMMA, 

GAMI, 

PI, 

1 IOPERXY, 

DIS2X, 

DIS2Y, 

DIS4X, 

DIS4Y, 

1 SMU, 

SMUIM, 

PHIDT, 

THETADT, 

METH, 

1 jacdt, 

DTRATE, 

nsuper. 

maxres(2), 

maxdq(2), 

1 RESID, 

RESIDMX, 

STRTIT, 

BCAIRF, 


1 CPUTIME, 

RESTART, 

STORE, 

IREAD, 

IPRINT, 

1 JTAIL1, 

JTAIL2, 

dswall, 

sobmax, 

CIRCUL, 

1 SHARP, 

CUSP, 

totime, 



1 numiter, 

i start, 

ISPEC, 

re, 

VISCOUS, 

1 IVIS, 

TURBULNT, VISXI, 

VISETA, 

VISCROSS, 

1 TRANSUP, 

TRANSLO, 

JTRANUP, 

JTRANLO, 

NPCP 


LOGICAL RESTART, STORE, TURBULNT, VISCOUS, PERIODIC, CIRCUL 
LOGICAL SHARP, CUSP, BCAIRF 
LOGICAL VISXI,VISETA,VISCROSS 
REAL*4 MXLNGTH 

COMMON /TEMPT/TINF,TWALL,WTRAT,TMN 


DIMENSION Q(JDIM,KDIM,4),TURMU(JDIM,KDIM),VORT(JDIM,KDIM) 
DIMENSION PRESS(JDIM,KDIM),XY(JDIM,KDIM,4),XYJ(JDIM,KDIM) 
DIMENSION X(JDIM,KDIM),Y(JDIM,KDIM) 


c 
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COMMONAVORKSP/SNOR(MAXJ),TMO(MAXJ),TMI(MAXJ),UU(MAXJ), 

* TAS(MAXJ),WORK(MAXJ,87) 

COMMON /FTURB/ FY(260,65),YMAXX(260) ( FMAXX(260),RYSM(260),YFMN 
COMMON/PLTDAT/ RESD(41000),CFPLT(MAXJ),CLPLT(4 1000), CDPLT(4 1000), 
&MXLNGTH(260,65),TMUI(260,65),TMUO(260,65), 

& DELTA(MAXJ) 

dimension FMUtmp( 2) 


C 

DATA VK,APLUS/0.4,26./ 
DATA Fl,F2/2000.,5.0/ 

C 

IF(TURBULNT) THEN 
IF(NUMITER.LT. 10) RETURN 
KEDGE = 0.75*KEND 
DO 40 J = JTAIL 1, JTAIL2 
C 

C DETERMINE VORTICITY TAS(K) 

C 

DO 5 K=KBEGIN,KUP 
TAS(K) = VORT(J,K) 
TURMU(J,K) = 0.0 
5 CONTINUE 
C 

C COMPUTE RA 
c 

do 15 k = KBEGIN,KBEGIN+1 
C2B = 198.6/TINF 
C2BP = C2B + 1. 
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RINV = 17Q(J,K,1) 

TT = GAMMA* PRESS(J,K)* RINV 
FMUtmp(K) = C2BP*TT**1.5/( C2B+TT) 

15 continue 
c 

K=KBEGIN 

WMU = .5*( FMUtmp(l) + FMUtmp(2)) 

TAU = 0. 10*ABS(VORT(J-2,K))+0.20*ABS(VORT(J-1,K» 

& + 0.40*ABS(VORT(J,K))+0.20*ABS(VORT(J+1,K)> 

& + 0.10*ABS(VORT(J+2,K)) 

RA = SQRT( RE*XYJ(J,K)*Q(J,K,1)*TAU/WMU) 

C WRITE(6,1000) TAU.RA 

C1000 FORMAT( IX, ’TAU = ‘,E10.4,2X,’RA = ‘.E10.4) 

C 

C CALCULATE NORMAL DISTANCE AND YSTAR 
C 

YSTAR = 1/RA 
DELTA(J) = YSTAR 
SNOR(l) = 0.0 
DO 10 K=KLOW,KUP 

SCIS = ABS(XY(J,K-1,3)*XY(J,K,3)+XY(J,K-1,4)*XY(J,K,4)) 
SCAL = 1.0/SQRT(SCIS) 

SNOR(K) = SNOR(K-l) + SCAL 
10 CONTINUE 
C 

C CALCULATE MIXING LENGTH 
C 

DO 30 K=KBEGIN,KEDGE 
YPLUS = SNOR(K)/YSTAR 
IF(YPLUS .LE. FI) THEN 

MXLNGTH(J,K) = VK*(Fl/F2)*YSTAR*(l.-(l.-YPLUS/Fl)**F2) 
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& *(1.-EXP(-YPLUS/APLUS)) 

ELSE 

MXLNGTH(J,K) = VK*(F1/F2)*YSTAR 
ENDIF 

TLSQR = MXLNGTH(J,K)*MXLNGTH(J,K) 

TURMU(J,K) = Q(J,K,1)*XYJ(J,K)*RE*TLSQR*ABS(V0RT(J,K)) 

30 CONTINUE 

40 CONTINUE 

C 

C USE CONSTANT T.E. LENGTH SCALE IN WAKE 
C 

DO 60 K=KLOW,KUP 
DO 50 J=1,JMAX 

IF(X(J,K) .GT. 1.0 .AND. Y(J,K) .LT. 0.0) JTELO=J 
IF(X(J,K) .LT. 1.0 .AND. Y(J,K) .GT. 0.0) JTEHI=J 
50 CONTINUE 

DO 60 J=1,JMAX 

IF(X(J,K) .GT. 1.0 .AND. Y(J,K) .LT. 0.0) THEN 
MXLNGTH(J.K) = MXLNGTH(JTELO.K) 

TLSQR = MXLNGTH(J,K)*MXLNGTH(J,K) 

TURMU(J.K) = Q(J,K,l)*XYJ(J,K)*RE*TLSQR*ABS(VORT(J,K)) 

ENDIF 

IF(X(J,K) .GT. 1.0 .AND. Y(J,K) .GT. 0.0) THEN 
MXLNGTH(J.K) = MXLNGTH(JTEHI.K) 

TLSQR = MXLNGTH(J,K)*MXLNGTH(J,K) 

TURMU(J.K) = Q(J,K,l)*XYJ(J,K)*RE*TLSQR + ABS(VORT(J,K)) 

ENDIF 

60 CONTINUE 

C 

C* ++ * PROVIDE FOR NONZERO LENGTH SCALE ACROSS CENTERLINE OF C-GRID 
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c 


JT1M1 = JTAIL1-1 
DO 70 J=1,JT1M1 
JMX2 = JMAX-J+1 

SLOPEML = (MXLNGTH(JMX2,5)-MXLNGTH(J,5))/(Y(JMX2,5)-Y(J,5)) 

DO 70 K=l,4 
C 

C**** CALCULATE MIXING LENGTH BELOW CENTERLINE 
C 

MXLNGTH(J.K) = SLOPEML*(Y(J,K)-Y(J,5))+MXLNGTH(J,5) 

TLSQR = MXLNGTH(J,K)*MXLNGTH(J,K) 

TURMU(J.K) = Q( J,K, 1 )*XYJ( J, K)*RE*TLSQR* ABS(V 0 RT( J,K)) 

C 

C**** CALCULATE MIXING LENGTH ABOVE CENTERLINE 
C 

MXLNGTH(JMX2,K) = MXLNGTH(JMX2,5)-SLOPEML*(Y(JMX2,5)-Y(JMX2,K)) 
TLSQR = MXLNGTH(JMX2,K)*MXLNGTH(JMX2,K) 

TURMU(JMX2,K)=Q(JMX2,K,l)*XYJ(JMX2,K)*RE*TLSQR*ABS(VORT(JMX2,K)) 
70 CONTINUE 


ifITRANSLO.ne.O.O)then 

c 

c zero turmu from JTRANLO to JTRANUP 


do 455 j = JTRANLO , JTRANUP 
do 455 k = KBEGIN.KEND 
turmu(j,k) = 0. 

455 continue 
c 

endif 
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c 


C IF NOT TURBULNT SET TURMU = 0.0 
C 

ELSE 

DO 800 K = KBEGIN.KEND 
DO 800 J = JBEGIN.JEND 
TURMU(J.K) = 0.0 
800 CONTINUE 
C 

ENDIF 

C 

RETURN 

END 


Appendix 3 


Coordinates of the Airfoil/Ice-Shape Geoemtries 

This appendix contains the x-y coordinate pairs that define the surface of the three 
airfoil/ice-shape geometries used in this work. The coordinates are listed as they would appear 
for an input file to the GRAPE grid generation code. All the x-coordinates are listed first 
with their corresponding y-values listed afterward. Both sets are listed in the same order(i.e. 
starting at the trailing edge and proceeding in a clockwise direction until reaching the trailing 
edge again.) All coordinates have been normalized by the chord length with the leading edge 
at (0.0, 0.0) and the trailing edge at (1.0, 0.0). Since the ice shapes grow in the negative x- 
direction, some coordinate locations will have x-values less than zero. 
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Coordinates for the NACA 63^- 415 airfoil and R7 ice shape 


X= 


1.0000000, 

0.9497200, 

0.8994100, 

0.6989400, 

0.6490700, 

0.5993000, 

0.4009500, 

0.3514800, 

0.3020000, 

0.1035300, 

0.0785300, 

0.0534000, 

0.0037000, 

-0.0013900, 

-0.0077790, 

-0.0081500, 

-0.0041700, 

0.0000000, 

0.0219800, 

0.0466000, 

0.0714700, 

0.2475000, 

0.2980000, 

0.3485200, 

0.5503900, 

0.6070000, 

0.6509300, 

0.8508500, 

0.9005900, 

0.9502800, 

Y= 



0.0000000, 

0.0033300, 

0.0018400, 

-0.0198900, 

-0.0266000, 

-0.0331100, 

-0.0524300, 

-0.0543900, 

-0.0547400, 

-0.0400900, 

-0.0356500, 

-0.0300000, 

-0.0181500, 

-0.0175900, 

-0.0168500, 

0.0000000, 

0.0063000, 

0.0115700, 

0.0296400, 

0.0426400, 

0.0526100, 

0.0894100, 

0.0936200, 

0.0955900, 

0.0829800, 

0.0759500, 

0.0678000, 

0.0288500, 

0.0188400, 

0.0093100, 


0.8491500, 

0.7989800, 

0.7489100, 

0.5496100, 

0.5000000, 

0.4504500, 

0.2525000, 

0.2029500, 

0.1533100, 

0.0280200, 

0.0150900, 

0.0099999, 

-0.0113000, 

-0.0131500, 

-0.0115700, 

0.0030000, 

0.0052500, 

0.0099100, 

0.0964700, 

0.1466900, 

0.1970500, 

0.3990500, 

0.4495500, 

0.5000000, 

0.7010600, 

0.7510900, 

0.8010200, 

1.0000000, 



-0.0019300, 

-0.0071600, 

-0.0132700, 

-0.0391800, 

-0.0445900, 

-0.0490900, 

-0.0536100, 

-0.0509500, 

-0.0465600, 

-0.0222000, 

-0.0164600, 

-0.0185200, 

-0.0151900, 

-0.0116700, 

-0.0060200, 

0.0128700, 

0.0158500, 

0.0207400, 

0.0607700, 

0.0734800, 

0.0827900, 

0.0952700, 

0.0928900, 

0.0887100, 

0.0587700, 

0.0490700, 

0.0390000, 

0.0000000 
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Coordinates for the NACA63 A -415 airfoil and G3 ice shape 


X= 


1.0000000, 

0.9497200, 

0.8994100, 

0.6989400, 

0.6490700, 

0.5993000, 

0.4009500, 

0.3514800, 

0.3020000, 

0.1035300, 

0.0785300, 

0.0534000, 

0.0065130, 

0.0046010, 

0.0022330, 

-0.0069620, 

-0.0081510, 

-0.0085260, 

-0.0049140, 

-0.0046560, 

-0.0048770, 

-0.0069890, 

-0.0077240, 

-0.0084540, 

-0.0150200, 

-0.0166600, 

-0.0184100, 

-0.0112500, 

-0.0070600, 

-0.0030580, 

0.0106100, 

0.0219800, 

0.0466000, 

0.1970500, 

0.2475000, 

0.2980000, 

0.5000000, 

0.5503900, 

0.6070000, 

0.8010200, 

0.8508500, 

0.9005900, 

Y= 



0.0000000, 

0.0033300, 

0.0018400, 

-0.0198900, 

-0.0266000, 

-0.0331100, 

-0.0524300, 

-0.0543900, 

-0.0547400, 

-0.0400900, 

-0.0356500, 

-0.0300000, 

-0.0139000, 

-0.0138300, 

-0.0135900, 

-0.0132400, 

-0.0126400, 

-0.0116800, 

-0.0090120, 

-0.0076940, 

-0.0040130, 

0.0061400, 

0.0068340, 

0.0070020, 

0.0079010, 

0.0088490, 

0.0105000, 

0.0128500, 

0.0128000, 

0.0129300, 

0.0195600, 

0.0296400, 

0.0426400, 

0.0827900, 

0.0894100, 

0.0936200, 

0.0887100, 

0.0829800, 

0.0747735, 

0.0390000, 

0.0288500, 

0.0188400, 


0.8491500, 

0.7989800, 

0.7489100, 

0.5496100, 

0.5000000, 

0.4504500, 

0.2525000, 

0.2029500, 

0.1533100, 

0.0280200, 

0.0116200, 

0.0088810, 

-0.0000456, 

-0.0028660, 

-0.0046860, 

-0.0079870, 

-0.0070810, 

-0.0055410, 

-0.0051820, 

-0.0056810, 

-0.0062490, 

-0.0094570, 

-0.0111000, 

-0.0132800, 

-0.0187900, 

-0.0178000, 

-0.0155300, 

0.0003970, 

0.0030300, 

0.0054740, 

0.0714700, 

0.0964700, 

0.1466900, 

0.3485200, 

0.3990500, 

0.4495500, 

0.6509300, 

0.7010600, 

0.7510900, 

0.9502800, 

1.0000000, 


-0.0019300, 

-0.0071600, 

-0.0132700, 

-0.0391800, 

-0.0445900, 

-0.0490900, 

-0.0536100, 

-0.0509500, 

-0.0465600, 

-0.0222000, 

-0.0148100, 

-0.0140500, 

-0.0132600, 

-0.0132900, 

-0.0133100, 

-0.0109700, 

-0.0106100, 

-0.0099830, 

-0.0010340, 

0.0029950, 

0.0050950, 

0.0071670, 

0.0072390, 

0.0073050, 

0.0120700, 

0.0130500, 

0.0131600, 

0.0132300, 

0.0137800, 

0.0149500, 

0.0526100, 

0.0607700, 

0.0734800, 

0.0955900, 

0.0952700, 

0.0928900, 

0.0678000, 

0.0587700, 

0.0490700, 

0.0093100, 

0.0000000 
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Coordinates for the NACA 0012 airfoil and G1 ice shape 




1.0000000, 

0.9913444, 

0.9809285, 

0.9681920, 

0.9523343, 

0.9321911, 

0.9060580, 

0.8714901, 

0.8252681, 

0.7642254, 

0.6882111, 

0.0043302, 

0.5252967, 

0.4601904, 

0.4103649, 

0.3700089, 

0.3300074, 

0.2887478, 

0.2051575, 

0.0888420, 

0.0851143, 

0.0737720, 

0.0695000, 

0.0657200, 

0.0619175, 

0.0581143, 

0.0543111, 

0.0505079, 

0.0467048, 

0.0390984, 

0.0352952, 

0.0314291, 

0.0309070, 

0.0276889, 

0.0238857, 

0.0233802, 

0.0228841, 

0.0224050, 

0.0219731, 

0.0216107, 

0.0213361, 

0.0211034, 

0.0188905, 

0.0166191, 

0.0143476, 

0.0120762, 

0.0098048, 

0.0075333, 

0.0052619, 

0.0029905, 

0.0007191, 

-0.0015524, 

-0.0038238, 

-0.0060952, 

-0.0083667, 

-0.0106381, 

-0.0129095, 

-0.0151810, 

-0.0174524, 

-0.0197238, 

-0.0219952, 

-0.0242667, 

-0.0265381, 

-0.0268877, 

-0.0269931, 

-0.0208443, 

-0.0264583, 

-0.0258790, 

-0.0251724, 

-0.0244191, 

-0.0213483, 

-0.0182775, 

-0.0152068, 

-0.0121361, 

-0.0090653, 

-0.0059946, 

-0.0029238, 

0.000 1409, 

0.0032177, 

0.0062884, 

0.0093592, 

0.0124299, 

0.0155007, 

0.0185714, 

0.0737220, 

0.0836666, 

0.0965504, 

0.1080990, 

0.1285863, 

0.1792650, 

0.2023058, 

0.3143436, 

0.4387820, 

0.5316864, 

0.5894307, 

0.6300145, 

0.6700144, 

0.7106869, 

0.7707552, 

0.8509442, 

0.9156263, 

0,9540105, 

0.9764334, 

0.9905078, 

1.0000000, 






X ~ 

0.0000000, 

-0.0012323, 

-0.0026976, 

-0.0044639, 

- 0.0066250, 

-0.0093109, 

-0.0127008, 

-0.0170272, 

-0.0225450, 

- 0.0293820, 

-0.0371858, 

-0.0448207, 

-0.0509522, 

-0.0550485, 

-0.0574780, 

-0.0589117, 

-0.0597804, 

-0.0000007, 

-0.0577608, 

-0.0449812, 

-0.0441132, 

-0.0417323, 

-0.0407000, 

-0.0410710, 

-0.0426421, 

-0.0436131, 

-0.0445841, 

- 0.0455552, 

-0.0465262, 

-0.0474972, 

-0.0484683, 

-0.0494393, 

-0.0504103, 

-0.0513813, 

-0.0523524, 

-0.0524220, 

-0.0523769, 

-0.0522204, 

-0.0519603, 

-0.0516100, 

-0.0511872, 

-0.0507137, 

-0.0490488, 

-0.0472992, 

-0.0454587, 

-0.0435218, 

-0.0414819, 

-0.0393315, 

-0.0370620, 

-0.0346632, 

-0.0321230, 

-0.0294269, 

-0.0265574, 

-0.0234920, 

-0.0202050, 

-0.0166590, 

-0.0128070, 

-0.0085832, 

-0.0038910, 

0.001 1202, 

0.0076136, 

0.0152347, 

0.0259245, 

0.0266049, 

0.0269931, 

0.0280895, 
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0.0287365, 

0.0292182, 

0.0294798, 

0.0280588, 

0.0275816, 

0.0271044, 

0.0251956, 

0.0247184, 

0.0242412, 

0.0417323, 

0.0439715, 

0.0462317, 

0.0576065, 

0.0599511, 

0.0561746, 

0.0389335, 

0.0349618, 

0.0286748, 

0.0033242, 

0.0013507, 

0.0000000 


0.0294914, 

0.0290133, 

0.0285361, 

0.0266272, 

0.0261500, 

0.0256728, 

0.0237640, 

0.0232867, 

0.0228095, 

0.0481440, 

0.0511603, 

0.0560180, 

0.0505007, 

0.0460672, 

0.0426025, 

0.0195168, 

0.0114718, 

0.0063985, 
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